This supplementary R Markdown document provides the full analytical workflow and visual outputs referenced in the main manuscript “The legacy of past heatwaves on ‘off-host’ parasite stages: Reduced infection risk and costs in the Daphnia-Pasteuria system”. All analyses were performed in R 4.4.2 (2024-10-31). The dataset and scripts are provided to ensure full reproducibility.

1. Data preparation and processing

This section outlines the steps taken to clean, transform and structure our dataset for subsequent analysis. Our dataset comprises various measurements of a sample of biological entities, including spore counts, infection statuses and life cycle metrics such as clutch size, offspring numbers and developmental stages. However, the raw dataset contained missing values that needed to be addressed prior to analysis.

1.1 Column descriptions

Below is a detailed explanation of each column in the dataset. This breakdown is intended to provide clarity on the variable names, units, and biological significance.

  • ID: A unique identifier for each sample, used to distinguish individual observations.

  • Temperature: The experimental temperature (in °C) at which each sample was maintained. The possible values are 20, 30 and 40 °C, which were used to examine the effects of temperature on biological parameters.

  • Isolate: Identifier for the specific parasite genotype used in the experiment (e.g. C1, C14, C18, C24). Each genotype represents a distinct biological strain or lineage of the parasite.

  • Clutch: The number of clutch events recorded for each sample, representing reproductive events or cycles.

  • Nb_offspring: The total number of offspring produced by each sample across clutch events.

  • FirstClutchAge: The age (in days) at which the first clutch was recorded.

  • LastClutchAge: The age (in days) at which the last clutch was recorded.

  • DeathAge: Age (in days) at death, providing a measure of lifespan for each sample.

  • Accidental_death: Indicates whether the death was accidental. Missing values suggest non-accidental or unspecified causes.

  • Volume_counted: The volume of the sample (in microlitres) counted under observation and used to standardise spore counts.

  • InfectedStatus: The infection status of each sample, with possible values including ‘infected’, ‘not_inf’ (not infected) or ‘lost’ for missing data.

  • Number of mature spores: The raw count of mature spores in the sample.

  • Number of immature spores (imm): The raw count of immature spores representing early-stage spore development.

  • Nb of white spores (whites): The number of white spores, which are mature but may appear white instead of red at the microscope.

  • cauliflowers: A categorical indicator (yes/no) of whether cauliflower-like structures were observed in the sample.

  • grapes: A categorical indicator (yes/no) of the presence of grape-like structures, which are another morphological feature observed in spore samples.

  • mature_load: The calculated load of mature spores, standardised based on ‘Volume_counted’. This represents the density of mature spores in the sample volume.

  • imm_load: The load of immature spores, calculated similarly to ‘mature_load’ for standardised comparison.

  • ratio: The ratio of mature spores to the total spore load (mature and immature). This metric provides an indication of spore maturation and infection severity.

  • Batch: A grouping variable created based on ID, representing experimental batches. This factor enables batch effects to be controlled in the analysis.

  • group: A combination of isolation and temperature conditions (e.g. C1840) used to categorise samples in specific experimental settings.

1.2 Key Data Processing Steps

  1. Variable Conversion: Several columns required type conversion to facilitate analysis. For example, ‘Volume_counted’ and ‘Nb of white spores’ were converted to numeric values. During this process, some ‘NA’ values were introduced due to incompatible data entries.

  2. Factorization: The categorical variables ‘cauliflowers’, ‘grapes’ and ‘group’ were converted into factor variables. Additionally, the ‘Isolate’ and ‘Temperature’ variables were explicitly factored to ensure clear categorical representation.

  3. Data Subsetting and Exclusion: Observations with specific conditions, such as an ID of 50, were excluded to maintain data integrity.

  4. Calculation of Derived Variables:

    • Spore Load Calculations: Two new variables were calculated: ‘mature_load’ and ‘imm_load’, representing the load of mature and immature spores respectively, standardised per unit volume.
    • Infection Ratio: The ratio variable was calculated as the proportion of mature spores relative to the total spore load to provide an indicator of infection severity.

  5. Grouping and Batching: Observations were assigned to batches for experimental grouping based on their ID. This enables batch effects to be tracked and controlled in downstream analysis.

##        ID         Temperature Isolate       Clutch        Nb_offspring   
##  Min.   :   1.0   20:129      C14: 92   Min.   : 1.000   Min.   :  0.00  
##  1st Qu.: 169.0   30:125      C1 : 87   1st Qu.: 1.000   1st Qu.:  0.00  
##  Median : 627.0   40:143      C18:116   Median : 5.000   Median : 10.00  
##  Mean   : 517.2               C24:102   Mean   : 5.666   Mean   : 27.48  
##  3rd Qu.: 903.0                         3rd Qu.:10.000   3rd Qu.: 57.25  
##  Max.   :1048.0                         Max.   :16.000   Max.   :111.00  
##                                         NA's   :101      NA's   :1       
##  FirstClutchAge  LastClutchAge      DeathAge     Accidental_death
##  Min.   : 1.00   Min.   :10.00   Min.   :15.00   Min.   :1       
##  1st Qu.:14.00   1st Qu.:21.00   1st Qu.:39.00   1st Qu.:1       
##  Median :18.00   Median :42.00   Median :52.00   Median :1       
##  Mean   :20.41   Mean   :42.49   Mean   :49.26   Mean   :1       
##  3rd Qu.:24.00   3rd Qu.:67.00   3rd Qu.:67.00   3rd Qu.:1       
##  Max.   :67.00   Max.   :67.00   Max.   :67.00   Max.   :1       
##  NA's   :102     NA's   :107                     NA's   :391     
##  Volume_counted   InfectedStatus     Number of mature spores
##  Min.   : 150.0   Length:397         Length:397             
##  1st Qu.: 150.0   Class :character   Class :character       
##  Median : 150.0   Mode  :character   Mode  :character       
##  Mean   : 274.9                                             
##  3rd Qu.: 450.0                                             
##  Max.   :1200.0                                             
##  NA's   :12                                                 
##  Number of imature spores Nb of white spores                cauliflowers
##  Length:397               Length:397         counted-after-freezer:  7  
##  Class :character         Class :character   no                   :270  
##  Mode  :character         Mode  :character   yes                  : 98  
##                                              NA's                 : 22  
##                                                                         
##                                                                         
##                                                                         
##                    grapes        whites            group      statusDeath  
##  counted-after-freezer:  7   Min.   :  0.000   C1840  : 43   Min.   :0.00  
##  no                   :266   1st Qu.:  0.000   C2420  : 39   1st Qu.:0.00  
##  yes                  :101   Median :  0.000   C1440  : 38   Median :1.00  
##  NA's                 : 23   Mean   :  6.753   C1830  : 38   Mean   :0.67  
##                              3rd Qu.:  1.000   C1820  : 35   3rd Qu.:1.00  
##                              Max.   :348.000   C2440  : 34   Max.   :1.00  
##                                                (Other):170                 
##  Number of immature spores  Never_repro         Batch          imm        
##  Length:397                Min.   :0.0000   11     : 12   Min.   :  0.00  
##  Class :character          1st Qu.:0.0000   6      : 11   1st Qu.:  0.00  
##  Mode  :character          Median :0.0000   28     : 11   Median :  0.00  
##                            Mean   :0.2544   4      : 10   Mean   : 16.58  
##                            3rd Qu.:1.0000   8      : 10   3rd Qu.: 18.00  
##                            Max.   :1.0000   14     : 10   Max.   :233.00  
##                                             (Other):333                   
##       mat          mature_load          imm_load           ratio       
##  Min.   :  0.00   Min.   :       0   Min.   :      0   Min.   :0.0000  
##  1st Qu.:  0.00   1st Qu.:       0   1st Qu.:      0   1st Qu.:0.7821  
##  Median :  0.00   Median :       0   Median :      0   Median :0.8734  
##  Mean   : 96.77   Mean   : 2357630   Mean   : 382870   Mean   :0.8179  
##  3rd Qu.:198.00   3rd Qu.: 4207500   3rd Qu.: 360000   3rd Qu.:0.9414  
##  Max.   :899.00   Max.   :20227500   Max.   :4230000   Max.   :1.0000  
##                   NA's   :12         NA's   :12        NA's   :250     
##    repro_rate               stade_max  
##  Min.   :0.0000   None           :235  
##  1st Qu.:0.0000   Cauliflower    : 14  
##  Median :0.2440   Grapes         :  1  
##  Mean   :0.4607   Immature       :  2  
##  3rd Qu.:0.9254   Mature         : 68  
##  Max.   :1.6567   2nd_Cauliflower: 77  
##  NA's   :1

1.3 Objectives and Methodology

  1. Objective: primary goal is to determine the effect of genotype (Isolate) and temperature conditions (Temperature) on each of the measured traits. The aim is to assess the individual and interactive effects of these variables while accounting for the random effects associated with batch variability.

  2. Model Selection:

    • Step 1: Random Effect Assessment – We fit two models to evaluate the importance of including batch as a random effect. We compare the model with the random effect against a fixed-only model using AICc to determine whether batch variability significantly affects infection probability.
    • Step 2: Fixed Effect Assessment – A series of models were fit with varying fixed effects (Isolate, Temperature, and their interaction), each including the batch as a random effect. The models were then evaluated using AICc to select the optimal fixed-effect structure.
    • Step 3: Final Model – Based on model selection criteria, the final model was fitted to estimate infection probability as a function of temperature while accounting for batch as a random effect. The performance of this model was assessed through residual simulation to confirm its suitability.

  3. Results and Interpretation: The final model outputs are summarised, enabling us to interpret the fixed effects of temperature on infection likelihood while controlling for batch variability.

2. Infectious Rate Analysis

This section examines the impact of temperature and the parasite’s genotype on infection rates. Specifically, we use generalised linear mixed models to examine how these factors, along with random batch effects, contribute to the probability of infection.

crushed$InfectedStatus = factor(crushed$InfectedStatus, levels=c("not_inf","infected"))

# Step 1: Evaluating the importance of random effect (Batch)
Inf5 = glmmTMB(data = crushed, InfectedStatus ~ Isolate * Temperature, family = "binomial", REML = TRUE)
Inf5R = glmmTMB(data = crushed, InfectedStatus ~ Isolate * Temperature + (1 | Batch), family = "binomial", REML = TRUE)

# Comparing models with and without random effect using AICc
AICc(Inf5, Inf5R) 
##       df     AICc
## Inf5  12 524.4920
## Inf5R 13 520.4055
# Step 2: Fixed effects selection
Inf0 = glmmTMB(data = crushed, InfectedStatus ~ 1 + (1 | Batch), family = "binomial", REML = FALSE)
Inf1 = glmmTMB(data = crushed, InfectedStatus ~ Isolate + (1 | Batch), family = "binomial", REML = FALSE)
Inf2 = glmmTMB(data = crushed, InfectedStatus ~ Temperature + (1 | Batch), family = "binomial", REML = FALSE)
Inf3 = glmmTMB(data = crushed, InfectedStatus ~ Isolate + Temperature + (1 | Batch), family = "binomial", REML = FALSE)
Inf4 = glmmTMB(data = crushed, InfectedStatus ~ Isolate * Temperature + (1 | Batch), family = "binomial", REML = FALSE)

# AICc comparison for fixed effects models
AICc(Inf0, Inf1, Inf2, Inf3, Inf4)
##      df     AICc
## Inf0  2 528.0016
## Inf1  5 531.1036
## Inf2  4 510.8766
## Inf3  7 514.1415
## Inf4 13 515.6191
# Model summaries for selected models
tab_model(Inf2, Inf3, Inf4, Inf1)
  InfectedStatus InfectedStatus InfectedStatus InfectedStatus
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
(Intercept) 0.76 0.58 – 0.98 0.037 0.63 0.39 – 1.02 0.062 0.63 0.39 – 1.03 0.067 0.60 0.38 – 0.95 0.029
Temperature [linear] 0.42 0.29 – 0.62 <0.001 0.42 0.29 – 0.62 <0.001 0.42 0.19 – 0.95 0.037
Temperature [quadratic] 1.03 0.71 – 1.50 0.866 1.05 0.72 – 1.53 0.801 0.63 0.29 – 1.40 0.262
Isolate [C1] 1.65 0.86 – 3.16 0.132 1.74 0.88 – 3.42 0.109 1.72 0.91 – 3.23 0.094
Isolate [C18] 1.23 0.67 – 2.26 0.497 1.21 0.66 – 2.22 0.545 1.28 0.71 – 2.30 0.415
Isolate [C24] 1.03 0.55 – 1.93 0.918 0.99 0.52 – 1.90 0.977 1.17 0.64 – 2.13 0.616
Isolate [C1] ×
Temperature [linear] 		0.68 0.21 – 2.22 0.519
Isolate [C18] ×
Temperature [linear] 		1.70 0.59 – 4.91 0.324
Isolate [C24] ×
Temperature [linear] 		0.62 0.20 – 1.91 0.400
Isolate [C1] ×
Temperature [quadratic] 		3.40 1.08 – 10.67 0.036
Isolate [C18] ×
Temperature [quadratic] 		2.13 0.74 – 6.13 0.162
Isolate [C24] ×
Temperature [quadratic] 		0.94 0.31 – 2.85 0.913
Random Effects
σ2 3.29 3.29 3.29 3.29
τ00 0.28 Batch 0.28 Batch 0.34 Batch 0.23 Batch
ICC 0.08 0.08 0.09 0.06
N 50 Batch 50 Batch 50 Batch 50 Batch
Observations 385 385 385 385
Marginal R2 / Conditional R2 0.066 / 0.139 0.076 / 0.149 0.117 / 0.201 0.010 / 0.073 
# Step 3: Final model selection based on AICc and interpretability
Inf_final = glmmTMB(data = crushed, InfectedStatus ~ Temperature + (1 | Batch), family = "binomial", REML = TRUE)

# Validating final model with residuals simulation
simulateResiduals(Inf_final, plot = TRUE)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2662472 0.337641 0.2329939 0.1066909 0.2742451 0.0246008 0.5366342 0.9108192 0.1563907 0.003480057 0.5733479 0.4877609 0.14512 0.3991937 0.304557 0.4133086 0.2767652 0.4752524 0.6629793 0.8902504 ...
# Summarizing the final model without intercept
tab_model(Inf_final, show.intercept = FALSE)
  InfectedStatus
Predictors Odds Ratios CI p
Temperature [linear] 0.44 0.30 – 0.63 <0.001
Temperature [quadratic] 1.03 0.71 – 1.49 0.872
Random Effects
σ2 3.29
τ00 Batch 0.30
ICC 0.08
N Batch 50
Observations 385
Marginal R2 / Conditional R2 0.061 / 0.140


Higher temperatures endured by spores significantly reduce infection rates, with the odds of infection decreasing by 57% per temperature unit increase (OR = 0.44, p < 0.001). No intermediate optimal temperature effect was observed, which suggests a consistent linear decrease. Although batch effects accounted for some variability, temperature remains the main driver.

Visualisation

3. Survival Analysis

This section explore how temperature and the parasite’s genotype contribute to the survival of infected individuals (also known as virulence) and uninfected individuals (potentially related to defense costs).

3.1 Infected individuals survival (= parasite virulence)

surv_obj = Surv(time = infected$DeathAge, event = infected$statusDeath)

fit_cox_5 = coxph(surv_obj ~ Isolate * Temperature, data = infected)
fit_cox_5R = coxme(surv_obj ~ Isolate * Temperature + (1|Batch), data = infected)

AICc(fit_cox_5, fit_cox_5R) # Random effect is significantly improving our results
##                  df     AICc
## fit_cox_5  11.00000 1329.365
## fit_cox_5R 28.43401 1326.367
fit_cox_1 = coxme(surv_obj ~ 1 + (1|Batch), data = infected)
fit_cox_2 = coxme(surv_obj ~ Temperature + (1|Batch), data = infected)
fit_cox_3 = coxme(surv_obj ~ Isolate + (1|Batch), data = infected)
fit_cox_4 = coxme(surv_obj ~ Isolate + Temperature + (1|Batch), data = infected)
fit_cox_5 = coxme(surv_obj ~ Isolate * Temperature + (1|Batch), data = infected)

AICc(fit_cox_1,fit_cox_2,fit_cox_3,fit_cox_4,fit_cox_5)
##                 df     AICc
## fit_cox_1 13.67456 1328.492
## fit_cox_2 16.10514 1331.733
## fit_cox_3 16.71636 1316.222
## fit_cox_4 19.33658 1319.417
## fit_cox_5 28.43401 1326.367
summary(fit_cox_3)
## Mixed effects coxme model
##  Formula: surv_obj ~ Isolate + (1 | Batch) 
##     Data: infected 
## 
##   events, n = 153, 167
## 
## Random effects:
##   group  variable        sd  variance
## 1 Batch Intercept 0.3965349 0.1572399
##                   Chisq    df         p   AIC    BIC
## Integrated loglik 21.73  4.00 2.270e-04 13.73   1.61
##  Penalized loglik 53.57 16.72 9.682e-06 20.14 -30.52
## 
## Fixed effects:
##               coef exp(coef) se(coef)     z      p
## IsolateC1   0.5158    1.6751   0.2618  1.97 0.0488
## IsolateC18 -0.5384    0.5837   0.2580 -2.09 0.0369
## IsolateC24 -0.1406    0.8688   0.2588 -0.54 0.5869
### Test assumptions of proportional hazards
(test_ph = cox.zph(fit_cox_3))
##         chisq df   p
## Isolate  6.15  3 0.1
## GLOBAL   6.15  3 0.1
plot(test_ph)

# Extract results from summary
summary_results = summary(fit_cox_3)

# Extract coefficients, standard errors, and p-values
coef = summary_results$coefficients[, "coef"]
se_coef = summary_results$coefficients[, "se(coef)"]
p_values = summary_results$coefficients[, "p"]

# Calculate HRs and confidence intervals
HR = exp(coef)
lower_CI = exp(coef - 1.96 * se_coef)
upper_CI = exp(coef + 1.96 * se_coef)

# Combine results into a data frame
data.frame(
  Predictors = rownames(summary_results$coefficients),
  HR = HR,
  Lower_95_CI = lower_CI,
  Upper_95_CI = upper_CI,
  p_value = p_values
)
##            Predictors        HR Lower_95_CI Upper_95_CI    p_value
## IsolateC1   IsolateC1 1.6750547   1.0026986   2.7982569 0.04880475
## IsolateC18 IsolateC18 0.5837033   0.3520170   0.9678783 0.03692993
## IsolateC24 IsolateC24 0.8688436   0.5232149   1.4427900 0.58690328

Visualisation

3.2 Exposed but non-infected Daphnia (potential defense costs)

surv_obj = Surv(time = not_infected$DeathAge, event = not_infected$statusDeath)

fit_cox_5 = coxph(surv_obj ~ Isolate * Temperature, data = not_infected)
fit_cox_5R = coxme(surv_obj ~ Isolate * Temperature + (1|Batch), data = not_infected)

AICc(fit_cox_5, fit_cox_5R) # Random effect
##                  df     AICc
## fit_cox_5  11.00000 1125.461
## fit_cox_5R 27.43641 1130.994
fit_cox_1 = coxph(surv_obj ~ 1 , data = not_infected)
fit_cox_2 = coxph(surv_obj ~ Temperature, data = not_infected)
fit_cox_3 = coxph(surv_obj ~ Isolate , data = not_infected)
fit_cox_4 = coxph(surv_obj ~ Isolate + Temperature , data = not_infected)
fit_cox_5 = coxph(surv_obj ~ Isolate * Temperature , data = not_infected)
AICc(fit_cox_1,fit_cox_2,fit_cox_3,fit_cox_4,fit_cox_5)
##           df     AICc
## fit_cox_1  0 1117.001
## fit_cox_2  2 1113.668
## fit_cox_3  3 1122.686
## fit_cox_4  5 1119.611
## fit_cox_5 11 1125.461
summary(fit_cox_2)
## Call:
## coxph(formula = surv_obj ~ Temperature, data = not_infected)
## 
##   n= 218, number of events= 110 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)   
## Temperature.L -0.4530    0.6357   0.1633 -2.773  0.00555 **
## Temperature.Q  0.1291    1.1379   0.1705  0.757  0.44882   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## Temperature.L    0.6357     1.5730    0.4616    0.8756
## Temperature.Q    1.1379     0.8789    0.8146    1.5894
## 
## Concordance= 0.577  (se = 0.026 )
## Likelihood ratio test= 7.45  on 2 df,   p=0.02
## Wald test            = 7.98  on 2 df,   p=0.02
## Score (logrank) test = 8.21  on 2 df,   p=0.02
anova(fit_cox_1, fit_cox_2)
## Analysis of Deviance Table
##  Cox model: response is  surv_obj
##  Model 1: ~ 1
##  Model 2: ~ Temperature
##    loglik  Chisq Df Pr(>|Chi|)  
## 1 -558.50                       
## 2 -554.78 7.4457  2    0.02416 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Test assumptions of proportional hazards
(test_ph =  cox.zph(fit_cox_2))
##             chisq df     p
## Temperature  6.48  2 0.039
## GLOBAL       6.48  2 0.039
plot(test_ph) ### Not so bad

### Weibull approach to compare results
fit_weibull1 = survreg(surv_obj ~ 1 , data = not_infected, dist = "weibull")
fit_weibull2 = survreg(surv_obj ~ Temperature, data = not_infected, dist = "weibull")
fit_weibull3 = survreg(surv_obj ~ Isolate, data = not_infected, dist = "weibull")
fit_weibull4 = survreg(surv_obj ~ Isolate + Temperature, data = not_infected, dist = "weibull")
fit_weibull5 = survreg(surv_obj ~ Isolate * Temperature, data = not_infected, dist = "weibull")
AICc(fit_weibull1,fit_weibull2,fit_weibull3,fit_weibull4,fit_weibull5)
##              df     AICc
## fit_weibull1  2 1212.177
## fit_weibull2  4 1207.967
## fit_weibull3  5 1217.837
## fit_weibull4  7 1213.874
## fit_weibull5 13 1219.663
# Standardi@zed residuals plot
residuals_std = residuals(fit_weibull2, type = "matrix")
qqnorm(residuals_std)
qqline(residuals_std, col = "red")

summary(fit_weibull2) # Same results, temperature is the main effect, with a significant decrease in mortality risk at higher temperatures (p = 0.003).
## 
## Call:
## survreg(formula = surv_obj ~ Temperature, data = not_infected, 
##     dist = "weibull")
##                 Value Std. Error     z       p
## (Intercept)    4.3576     0.0589 74.01 < 2e-16
## Temperature.L  0.2701     0.0931  2.90  0.0037
## Temperature.Q -0.0777     0.0960 -0.81  0.4182
## Log(scale)    -0.5755     0.0864 -6.66 2.7e-11
## 
## Scale= 0.562 
## 
## Weibull distribution
## Loglik(model)= -599.9   Loglik(intercept only)= -604.1
##  Chisq= 8.34 on 2 degrees of freedom, p= 0.015 
## Number of Newton-Raphson Iterations: 5 
## n= 218
# So I decided to keep the cox model despite a slight violation of the proportional hazard assumption, as it is more interpretable and the results are consistent with the Weibull model.

# Extract results from summary
summary_results = summary(fit_cox_2)

# Extract coefficients, standard errors, and p-values
coef = summary_results$coefficients[, "coef"]
se_coef = summary_results$coefficients[, "se(coef)"]
p_values = summary_results$coefficients[, "Pr(>|z|)"]

# Calculate HRs and confidence intervals
HR = exp(coef)
lower_CI = exp(coef - 1.96 * se_coef)
upper_CI = exp(coef + 1.96 * se_coef)

# Combine results into a data frame
data.frame(
  Predictors = rownames(summary_results$coefficients),
  HR = HR,
  Lower_95_CI = lower_CI,
  Upper_95_CI = upper_CI,
  p_value = p_values
)
##                  Predictors        HR Lower_95_CI Upper_95_CI     p_value
## Temperature.L Temperature.L 0.6357221   0.4615649   0.8755921 0.005547673
## Temperature.Q Temperature.Q 1.1378501   0.8146002   1.5893723 0.448821505


The model indicates that temperature significantly decreases the hazard rate, with a 37.1% reduction in mortality risk (HR = 0.63, p = 0.006), suggesting potential temperature-mediated effects on host survival.

Visualisation

3.3 Exposed uninfected Daphnia (including unexposed husbandry control)

##   C1  C14  C18  C24 CTRL 
##   45   58   67   60   27

We start first by running an analysis including the husbandry control along with the other genotype.

Warning! It is important to note that under these conditions it is not possible to test the interaction between heatwave temperature and the absence of spores, as our husbandry controls were infected with crushed Daphnia bodies that did not experience any heatwaves.

surv_obj = Surv(time = dataC$DeathAge, event = dataC$statusDeath)

fit_cox_5 = coxph(surv_obj ~  Temperature + Isolate, data = dataC)
fit_cox_5R = coxme(surv_obj ~  Temperature + Isolate + (1|Batch), data = dataC) 
# Cannot evaluate the interaction as the unexposed individuals do not have a heatwave condition (no spores)

AICc(fit_cox_5, fit_cox_5R) # No random effect
##                  df     AICc
## fit_cox_5   6.00000 1376.118
## fit_cox_5R 12.21118 1376.970
fit_cox_1 = coxph(surv_obj ~ 1, data = dataC)
fit_cox_2 = coxph(surv_obj ~ Temperature, data = dataC)
fit_cox_3 = coxph(surv_obj ~ Isolate, data = dataC)
fit_cox_4 = coxph(surv_obj ~ Isolate + Temperature, data = dataC)

AICc(fit_cox_1,fit_cox_2,fit_cox_3,fit_cox_4)
##           df     AICc
## fit_cox_1  0 1376.198
## fit_cox_2  2 1368.589
## fit_cox_3  4 1379.257
## fit_cox_4  6 1376.118
### Test assumptions of proportional hazards
test_ph = cox.zph(fit_cox_2)
test_ph
##             chisq df    p
## Temperature  4.38  2 0.11
## GLOBAL       4.38  2 0.11
plot(test_ph) # Not the best. fit

### Weibull approach to compare results
fit_weibull1 = survreg(surv_obj ~ 1 , data = dataC, dist = "weibull")
fit_weibull2 = survreg(surv_obj ~ Temperature, data = dataC, dist = "weibull")
fit_weibull3 = survreg(surv_obj ~ Isolate, data = dataC, dist = "weibull")
fit_weibull4 = survreg(surv_obj ~ Isolate + Temperature, data = dataC, dist = "weibull")
fit_weibull5 = survreg(surv_obj ~ Isolate * Temperature, data = dataC, dist = "weibull")
AICc(fit_weibull1,fit_weibull2,fit_weibull3,fit_weibull4,fit_weibull5)
##              df     AICc
## fit_weibull1  2 1438.865
## fit_weibull2  4 1429.968
## fit_weibull3  6 1441.145
## fit_weibull4  8 1437.353
## fit_weibull5 16 1447.669
# Standardi@zed residuals plot
residuals_std = residuals(fit_weibull2, type = "matrix")
qqnorm(residuals_std)
qqline(residuals_std, col = "red")

summary(fit_weibull2) # not much much better, but sit's also because the factor are not equilibrated as the design is not meant to test those hypothesese.
## 
## Call:
## survreg(formula = surv_obj ~ Temperature, data = dataC, dist = "weibull")
##                 Value Std. Error     z       p
## (Intercept)    4.3525     0.0523 83.27 < 2e-16
## Temperature.L  0.2662     0.0783  3.40 0.00067
## Temperature.Q -0.1057     0.0888 -1.19 0.23395
## Log(scale)    -0.6169     0.0782 -7.89 3.1e-15
## 
## Scale= 0.54 
## 
## Weibull distribution
## Loglik(model)= -710.9   Loglik(intercept only)= -717.4
##  Chisq= 13.01 on 2 degrees of freedom, p= 0.0015 
## Number of Newton-Raphson Iterations: 5 
## n= 257
summary(fit_cox_2)
## Call:
## coxph(formula = surv_obj ~ Temperature, data = dataC)
## 
##   n= 257, number of events= 132 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)   
## Temperature.L -0.4683    0.6261   0.1435 -3.264   0.0011 **
## Temperature.Q  0.1847    1.2028   0.1644  1.123   0.2614   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## Temperature.L    0.6261     1.5973    0.4726    0.8294
## Temperature.Q    1.2028     0.8314    0.8715    1.6601
## 
## Concordance= 0.582  (se = 0.023 )
## Likelihood ratio test= 11.7  on 2 df,   p=0.003
## Wald test            = 12.24  on 2 df,   p=0.002
## Score (logrank) test = 12.64  on 2 df,   p=0.002

Secondly, we analysed the survival of uninfected Daphnias only at 20°C, as this is a condition were the spores endured no heatwave, like our husbandry controls.

##                 df     AICc
## fit_cox_5  4.00000 433.5022
## fit_cox_5R 4.01947 433.5129
##           df     AICc
## fit_cox_1  0 429.8015
## fit_cox_2  4 433.5022
## Call:
## coxph(formula = surv_obj ~ Isolate, data = dataC20)
## 
##   n= 84, number of events= 54 
## 
##               coef exp(coef) se(coef)      z Pr(>|z|)  
## IsolateC1   0.7564    2.1306   0.4438  1.704   0.0883 .
## IsolateC14 -0.1502    0.8605   0.4052 -0.371   0.7108  
## IsolateC18 -0.2094    0.8111   0.3910 -0.536   0.5923  
## IsolateC24 -0.4021    0.6689   0.4058 -0.991   0.3218  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## IsolateC1     2.1306     0.4694    0.8928     5.084
## IsolateC14    0.8605     1.1621    0.3889     1.904
## IsolateC18    0.8111     1.2329    0.3769     1.745
## IsolateC24    0.6689     1.4949    0.3020     1.482
## 
## Concordance= 0.571  (se = 0.044 )
## Likelihood ratio test= 5.12  on 4 df,   p=0.3
## Wald test            = 5.98  on 4 df,   p=0.2
## Score (logrank) test = 6.39  on 4 df,   p=0.2
##         chisq df    p
## Isolate  4.19  4 0.38
## GLOBAL   4.19  4 0.38


We can see that there is no significant difference between the model comparing the survival of each spore lineage (including the group without spores) and the null model. So it seem than at 20°C, in stable conditions, the exposure to dead Daphnias or to spores from different genotype influence similarly the survival of non infected individuals.

Visualisation


In conclusion, when we included the husbandry control in our analysis, comparing the additive effects of heatwave temperature and different spores genotypes (including non spores at all), the main effect explaining survival variation is still temperature.

4. Host reproduction

In this section we explored how the exposure and the infection by heat-stressed parasites affected the host reproduction. First we evaluated the number of individuals fully castrated in exposed uninfected and in infected Daphnias. This divisions into two sub analysis helped us to find more adapted models to describe the host reproductive patern.

4.1 Full castration

One of the consequences of Pasteuria ramosa infection in Daphnia is the reduction in allocation to reproduction, which can also potentially be due to the energetic cost of immune response for uninfected hosts. Here we focused on the proportion of non-reproducing individuals.

4.1.1 Infected individuals

FC1 = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge + (1|Batch), family = "binomial", REML=T)

AICc(FC1, FC1R)
##      df     AICc
## FC1   7 239.8814
## FC1R  8 238.7983
FC1 = glmmTMB(data=infected, Never_repro ~ Isolate * Temperature + DeathAge, family = "binomial", 
               na.action = na.fail)
model_selection = dredge(FC1)

# see best models
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge, 
##     data = infected, family = "binomial", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df   logLik
## 16    0.79900          + -0.02228        +        +            + 13  -98.284
## 15   -0.19180          +                 +        +            + 12  -99.622
## 4     0.98420          + -0.01982        +                        5 -108.183
## 3     0.11780          +                 +                        4 -109.410
## 2     0.97250          + -0.02725                                 2 -111.825
## 8     0.94670          + -0.01938        +        +               7 -107.700
## 7     0.09487          +                 +        +               6 -108.859
## 1    -0.25280          +                                          1 -114.432
## 6     0.91880          + -0.02678                 +               4 -111.381
## 5    -0.28880          +                          +               3 -113.875
##     AICc delta weight
## 16 224.9  0.00  0.321
## 15 225.3  0.32  0.273
## 4  226.7  1.79  0.131
## 3  227.1  2.12  0.111
## 2  227.7  2.78  0.080
## 8  230.1  5.16  0.024
## 7  230.2  5.30  0.023
## 1  230.9  5.94  0.016
## 6  231.0  6.06  0.016
## 5  233.9  8.95  0.004
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Never_repro Never_repro Never_repro
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
DeathAge 0.98 0.95 – 1.00 0.104 0.98 0.96 – 1.01 0.120
Isolate [C1] 2.27 0.66 – 7.83 0.193 2.46 0.72 – 8.38 0.151 1.23 0.49 – 3.10 0.658
Isolate [C18] 0.55 0.18 – 1.69 0.298 0.50 0.16 – 1.49 0.212 0.44 0.18 – 1.11 0.083
Isolate [C24] 0.65 0.20 – 2.09 0.472 0.66 0.21 – 2.08 0.474 0.50 0.20 – 1.28 0.150
Temperature [linear] 0.17 0.03 – 0.96 0.044 0.14 0.03 – 0.77 0.024
Temperature [quadratic] 0.28 0.07 – 1.12 0.073 0.29 0.07 – 1.16 0.080
Isolate [C1] ×
Temperature [linear] 		32.02 3.12 – 328.91 0.004 36.19 3.56 – 367.88 0.002
Isolate [C18] ×
Temperature [linear] 		2.04 0.26 – 16.14 0.497 2.80 0.37 – 21.15 0.319
Isolate [C24] ×
Temperature [linear] 		5.08 0.56 – 45.81 0.147 6.46 0.74 – 56.47 0.092
Isolate [C1] ×
Temperature [quadratic] 		10.98 1.60 – 75.39 0.015 10.14 1.50 – 68.66 0.018
Isolate [C18] ×
Temperature [quadratic] 		3.86 0.64 – 23.40 0.142 3.51 0.59 – 20.92 0.168
Isolate [C24] ×
Temperature [quadratic] 		3.68 0.59 – 23.04 0.163 3.31 0.54 – 20.28 0.196
Observations 167 167 167
R2 Tjur 0.178 0.163 0.074 
best_model = get.models(model_selection, 2)[[1]]
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
tab_model(best_model, show.intercept = F)
  Never_repro
Predictors Odds Ratios CI p
Isolate [C1] 2.46 0.72 – 8.38 0.151
Isolate [C18] 0.50 0.16 – 1.49 0.212
Isolate [C24] 0.66 0.21 – 2.08 0.474
Temperature [linear] 0.14 0.03 – 0.77 0.024
Temperature [quadratic] 0.29 0.07 – 1.16 0.080
Isolate [C1] ×
Temperature [linear] 		36.19 3.56 – 367.88 0.002
Isolate [C18] ×
Temperature [linear] 		2.80 0.37 – 21.15 0.319
Isolate [C24] ×
Temperature [linear] 		6.46 0.74 – 56.47 0.092
Isolate [C1] ×
Temperature [quadratic] 		10.14 1.50 – 68.66 0.018
Isolate [C18] ×
Temperature [quadratic] 		3.51 0.59 – 20.92 0.168
Isolate [C24] ×
Temperature [quadratic] 		3.31 0.54 – 20.28 0.196
Observations 167
R2 Tjur 0.163
Visualisation
## `summarise()` has grouped output by 'Temperature'. You can override using the
## `.groups` argument.


A heatwave experienced by the parasite appears to relax its castrating effect on the host for most genotypes, with a strong opposit effect for C1 genotyp, with a lower probability of total sterility in genotype C1, which shows a dramatically increased risk of host sterility after exposure to 40°C. This suggests an interaction between genotype and heatwave experienced by the parasite.

4.1.2 Uninfected individuals

options(na.action = "na.fail")

FC1 = glmmTMB(data=not_infected, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=not_infected, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)

AICc(FC1, FC1R)
##      df     AICc
## FC1   7 118.8705
## FC1R  8 106.9146
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)

# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge + 
##     (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0, 
##     dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df  logLik
## 2       1.672          +  -0.1281                                 3 -48.375
## 6       1.459          +  -0.1208                 +               5 -46.911
## 4       2.324          +  -0.1529        +                        6 -46.318
## 8       2.030          +  -0.1347        +        +               8 -45.341
## 16     13.600          +  -0.8717        +        +            + 14 -40.057
## 5      -2.909          +                          +               4 -67.021
## 7      -2.961          +                 +        +               7 -65.982
## 1      -2.830          +                                          2 -72.618
## 15     -3.165          +                 +        +            + 13 -61.318
## 3      -2.927          +                 +                        5 -71.734
##     AICc delta weight
## 2  102.9  0.00  0.499
## 6  104.1  1.24  0.268
## 4  105.0  2.17  0.168
## 8  107.4  4.51  0.052
## 16 110.2  7.32  0.013
## 5  142.2 39.37  0.000
## 7  146.5 43.63  0.000
## 1  149.3 46.43  0.000
## 15 150.4 47.56  0.000
## 3  153.8 50.89  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | Batch)
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Never_repro Never_repro
Predictors Odds Ratios CI p Odds Ratios CI p
DeathAge 0.88 0.82 – 0.94 <0.001 0.89 0.83 – 0.95 0.001
Temperature [linear] 0.36 0.10 – 1.30 0.120
Temperature [quadratic] 0.61 0.17 – 2.17 0.442
Random Effects
σ2 3.29 3.29
τ00 6.14 Batch 5.52 Batch
ICC 0.65 0.63
N 50 Batch 50 Batch
Observations 218 218
Marginal R2 / Conditional R2 0.385 / 0.786 0.421 / 0.784 
# best model
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
##  Family: binomial  ( logit )
## Formula:          Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.8    112.9    -48.4     96.8      215 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Batch  (Intercept) 6.137    2.477   
## Number of obs: 218, groups:  Batch, 50
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.67209    0.98865   1.691  0.09078 .  
## DeathAge    -0.12808    0.03582  -3.575  0.00035 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
  Never_repro
Predictors Odds Ratios CI p
(Intercept) 5.32 0.77 – 36.96 0.091
DeathAge 0.88 0.82 – 0.94 <0.001
Random Effects
σ2 3.29
τ00 Batch 6.14
ICC 0.65
N Batch 50
Observations 218
Marginal R2 / Conditional R2 0.385 / 0.786


Only the age at death seems to reprduce at least once in their life time.

Visualisation
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.

4.1.3 Uninfected individuals (including unexposed husbandry control)

options(na.action = "na.fail")

dataC_ninf = subset(dataC, InfectedStatus!="infected")
summary(dataC_ninf)
##        ID       Temperature Isolate       Clutch        Nb_offspring   
##  Min.   :   2   20: 84      C1  :45   Min.   : 1.000   Min.   :  0.00  
##  1st Qu.: 218   30: 73      C14 :58   1st Qu.: 3.750   1st Qu.: 10.00  
##  Median : 668   40:100      C18 :67   Median : 8.000   Median : 47.00  
##  Mean   : 588               C24 :60   Mean   : 7.368   Mean   : 43.02  
##  3rd Qu.: 969               CTRL:27   3rd Qu.:11.000   3rd Qu.: 68.00  
##  Max.   :1100                         Max.   :16.000   Max.   :111.00  
##                                       NA's   :29                       
##  FirstClutchAge  LastClutchAge      DeathAge     Accidental_death
##  Min.   : 1.00   Min.   :10.00   Min.   :15.00   Min.   :1       
##  1st Qu.:14.00   1st Qu.:39.00   1st Qu.:31.00   1st Qu.:1       
##  Median :18.00   Median :59.00   Median :62.00   Median :1       
##  Mean   :21.77   Mean   :50.47   Mean   :51.47   Mean   :1       
##  3rd Qu.:27.00   3rd Qu.:67.00   3rd Qu.:67.00   3rd Qu.:1       
##  Max.   :67.00   Max.   :67.00   Max.   :67.00   Max.   :1       
##  NA's   :30      NA's   :32                      NA's   :250     
##  Volume_counted InfectedStatus     Number of mature spores
##  Min.   :150    Length:257         Length:257             
##  1st Qu.:150    Class :character   Class :character       
##  Median :150    Mode  :character   Mode  :character       
##  Mean   :152                                              
##  3rd Qu.:150                                              
##  Max.   :200                                              
##  NA's   :12                                               
##  Number of imature spores Nb of white spores cauliflowers  grapes   
##  Length:257               Length:257         no  :236     no  :236  
##  Class :character         Class :character   NA's: 21     NA's: 21  
##  Mode  :character         Mode  :character                          
##                                                                     
##                                                                     
##                                                                     
##                                                                     
##      whites      group      statusDeath     Number of immature spores
##  Min.   :0   C1440  : 29   Min.   :0.0000   Length:257               
##  1st Qu.:0   C2440  : 28   1st Qu.:0.0000   Class :character         
##  Median :0   NACTRL : 27   Median :1.0000   Mode  :character         
##  Mean   :0   C1840  : 26   Mean   :0.5136                            
##  3rd Qu.:0   C1830  : 24   3rd Qu.:1.0000                            
##  Max.   :0   C130   : 19   Max.   :1.0000                            
##              (Other):104                                             
##   Never_repro         Batch          imm         mat     mature_load
##  Min.   :0.0000   2      :  8   Min.   :0   Min.   :0   Min.   :0   
##  1st Qu.:0.0000   17     :  8   1st Qu.:0   1st Qu.:0   1st Qu.:0   
##  Median :0.0000   28     :  8   Median :0   Median :0   Median :0   
##  Mean   :0.1128   31     :  8   Mean   :0   Mean   :0   Mean   :0   
##  3rd Qu.:0.0000   3      :  7   3rd Qu.:0   3rd Qu.:0   3rd Qu.:0   
##  Max.   :1.0000   6      :  7   Max.   :0   Max.   :0   Max.   :0   
##                   (Other):211                           NA's   :12  
##     imm_load      ratio       repro_rate     stade_max 
##  Min.   :0    Min.   : NA   Min.   :0.0000   None:257  
##  1st Qu.:0    1st Qu.: NA   1st Qu.:0.3095             
##  Median :0    Median : NA   Median :0.8358             
##  Mean   :0    Mean   :NaN   Mean   :0.7206             
##  3rd Qu.:0    3rd Qu.: NA   3rd Qu.:1.0597             
##  Max.   :0    Max.   : NA   Max.   :1.6567             
##  NA's   :12   NA's   :257
FC1 = glmmTMB(data=dataC_ninf, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=dataC_ninf, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)

AICc(FC1, FC1R)
##      df     AICc
## FC1   8 134.3687
## FC1R  9 125.3196
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)

# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge + 
##     (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0, 
##     dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df  logLik
## 2       1.672          +  -0.1281                                 3 -48.375
## 6       1.459          +  -0.1208                 +               5 -46.911
## 4       2.324          +  -0.1529        +                        6 -46.318
## 8       2.030          +  -0.1347        +        +               8 -45.341
## 16     13.600          +  -0.8717        +        +            + 14 -40.057
## 5      -2.909          +                          +               4 -67.021
## 7      -2.961          +                 +        +               7 -65.982
## 1      -2.830          +                                          2 -72.618
## 15     -3.165          +                 +        +            + 13 -61.318
## 3      -2.927          +                 +                        5 -71.734
##     AICc delta weight
## 2  102.9  0.00  0.499
## 6  104.1  1.24  0.268
## 4  105.0  2.17  0.168
## 8  107.4  4.51  0.052
## 16 110.2  7.32  0.013
## 5  142.2 39.37  0.000
## 7  146.5 43.63  0.000
## 1  149.3 46.43  0.000
## 15 150.4 47.56  0.000
## 3  153.8 50.89  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Never_repro Never_repro
Predictors Odds Ratios CI p Odds Ratios CI p
DeathAge 0.88 0.82 – 0.94 <0.001 0.89 0.83 – 0.95 0.001
Temperature [linear] 0.36 0.10 – 1.30 0.120
Temperature [quadratic] 0.61 0.17 – 2.17 0.442
Random Effects
σ2 3.29 3.29
τ00 6.14 Batch 5.52 Batch
ICC 0.65 0.63
N 50 Batch 50 Batch
Observations 218 218
Marginal R2 / Conditional R2 0.385 / 0.786 0.421 / 0.784 
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
##  Family: binomial  ( logit )
## Formula:          Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.8    112.9    -48.4     96.8      215 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Batch  (Intercept) 6.137    2.477   
## Number of obs: 218, groups:  Batch, 50
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.67209    0.98865   1.691  0.09078 .  
## DeathAge    -0.12808    0.03582  -3.575  0.00035 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
  Never_repro
Predictors Odds Ratios CI p
(Intercept) 5.32 0.77 – 36.96 0.091
DeathAge 0.88 0.82 – 0.94 <0.001
Random Effects
σ2 3.29
τ00 Batch 6.14
ICC 0.65
N Batch 50
Observations 218
Marginal R2 / Conditional R2 0.385 / 0.786 
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.9088821 0.9008778 0.9317442 0.7483929 0.241452 0.2361889 0.7475164 0.4387352 0.390151 0.5425091 0.609537 0.7492443 0.6104253 0.7648734 0.5350931 0.2879339 0.134203 0.1612245 0.3157523 0.1570396 ...


Only the age at death seems to reprduce at least once in their life time, the husbandry control inclusion does not change this.

Visualisation
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.

4.2 Number of individuals produced

In this section we were interested in understanding the role of parasite heat-stress and genotype on the number of offspring produced. Indeed, one of the consequences of infection in Daphnia is the reduction in allocation to parasite reproduction, (or a potential cost of immune response for uninfected individuals), even if the host is not fully castrated.

4.2.1 Infected individuals

repro_inf = subset(data, data$Clutch>0 & InfectedStatus=="infected" & Temperature!="NA")
repro_inf$Temperature = factor(repro_inf$Temperature, ordered = T)

FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "poisson", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=T)

simulateResiduals(FC1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.8741857 0.7996104 0.4441007 0.4213267 0.3856629 0.2460142 0.4371622 0.9950848 0.5756341 0.3335869 0.1840254 0.6387504 0.2376711 0.8736609 0.9228149 0.7387757 0.2361033 0.4922073 0.4105355 0.1956212 ...
FC1R = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "nbinom2", REML=T)
AICc(FC1, FC1R)# no batch effect
##      df     AICc
## FC1   8 601.8980
## FC1R  9 601.5721
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge  , 
                data = repro_inf, family = "nbinom2", na.action = na.fail)
model_selection = dredge(FC1)

model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge, 
##     data = repro_inf, family = "nbinom2", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df   logLik
## 4      0.9519          +  0.03153        +                        6 -283.958
## 2      0.2723          +  0.03680                                 3 -289.038
## 8      0.9455          +  0.03176        +        +               8 -283.902
## 16     0.7038          +  0.02896        +        +            + 14 -277.332
## 6      0.3050          +  0.03628                 +               5 -288.741
## 3      2.7120          +                 +                        5 -291.710
## 7      2.7390          +                 +        +               7 -291.682
## 15     2.2290          +                 +        +            + 13 -284.113
## 1      2.1350          +                                          2 -300.042
## 5      2.1400          +                          +               4 -299.296
##     AICc delta weight
## 4  580.9  0.00  0.750
## 2  584.3  3.46  0.133
## 8  585.5  4.62  0.075
## 16 588.0  7.10  0.022
## 6  588.2  7.28  0.020
## 3  594.1 13.22  0.001
## 7  598.7 17.79  0.000
## 15 598.8 17.89  0.000
## 1  604.2 23.34  0.000
## 5  607.0 26.16  0.000
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Nb_offspring
Predictors Incidence Rate Ratios CI p
DeathAge 1.03 1.02 – 1.05 <0.001
Isolate [C1] 0.44 0.22 – 0.89 0.021
Isolate [C18] 0.48 0.28 – 0.82 0.008
Isolate [C24] 0.80 0.45 – 1.39 0.424
Observations 94 
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T) 

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.8445634 0.7757172 0.3188654 0.4678265 0.4046863 0.2574613 0.3273075 0.988557 0.6029993 0.2770807 0.2222652 0.5937299 0.200338 0.8887268 0.8925955 0.6897923 0.2627872 0.457771 0.3542096 0.1563993 ...
tab_model(best_model, show.intercept = F)
  Nb_offspring
Predictors Incidence Rate Ratios CI p
DeathAge 1.03 1.02 – 1.05 <0.001
Isolate [C1] 0.44 0.22 – 0.89 0.021
Isolate [C18] 0.48 0.28 – 0.82 0.008
Isolate [C24] 0.80 0.45 – 1.39 0.424
Observations 94 
summary(best_model)
##  Family: nbinom2  ( log )
## Formula:          Nb_offspring ~ DeathAge + Isolate + 1
## Data: repro_inf
## 
##      AIC      BIC   logLik deviance df.resid 
##    579.9    595.2   -284.0    567.9       88 
## 
## 
## Dispersion parameter for nbinom2 family ():  1.4 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.951947   0.468983   2.030  0.04238 *  
## DeathAge     0.031525   0.007622   4.136 3.54e-05 ***
## IsolateC1   -0.812114   0.352549  -2.304  0.02125 *  
## IsolateC18  -0.738421   0.276490  -2.671  0.00757 ** 
## IsolateC24  -0.228773   0.286424  -0.799  0.42445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Visualisation

4.2.2 Uninfected individuals

repro_nf = subset(data, data$Clutch>0 & InfectedStatus=="not_inf" & Temperature!="NA")
repro_nf$Temperature = factor(repro_nf$Temperature, ordered = T)

FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "gaussian", REML=T)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature+ DeathAge + (1|Batch), family = "poisson")
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_nf, log(Nb_offspring) ~ Isolate * Temperature + DeathAge + (1|Batch), family = "gaussian", REML=T)

FC1R = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge  + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nf, family = "gaussian", REML=F)

AICc(FC1, FC1R)# Batch effect NOT important! But significantly improved the residuals
##      df     AICc
## FC1   8 1659.382
## FC1R  9 1660.230
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge + (1|Batch), data = repro_nf, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1)

model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge + 
##     (1 | Batch), data = repro_nf, family = "gaussian", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df   logLik
## 2      -22.63          +    1.306                                 4 -823.741
## 4      -20.37          +    1.313        +                        7 -821.201
## 6      -22.45          +    1.294                 +               6 -823.160
## 8      -20.41          +    1.303        +        +               9 -820.623
## 5       48.57          +                          +               5 -913.919
## 1       50.26          +                                          3 -916.804
## 7       51.89          +                 +        +               8 -913.084
## 3       53.77          +                 +                        6 -915.627
## 15      52.47          +                 +        +            + 14 -912.123
## 16     -20.73          +    1.307        +        +            + 15         
##      AICc  delta
## 2  1655.7   0.00
## 4  1657.0   1.31
## 6  1658.8   3.08
## 8  1660.2   4.54
## 5  1838.2 182.46
## 1  1839.7 184.04
## 7  1843.0 187.26
## 3  1843.7 188.01
## 15 1854.6 198.91
## 16              
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
selected_models
## $`2`
## Formula:          Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
##       AIC       BIC    logLik  df.resid 
## 1655.4818 1668.5326 -823.7409       189 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups   Name        Std.Dev.
##  Batch    (Intercept)  4.599  
##  Residual             16.721  
## 
## Number of obs: 193 / Conditional model: Batch, 50
## 
## Dispersion estimate for gaussian family (sigma^2):  280 
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)     DeathAge  
##     -22.631        1.306  
## 
## $`4`
## Formula:          Nb_offspring ~ DeathAge + Isolate + (1 | Batch)
## Data: repro_nf
##       AIC       BIC    logLik  df.resid 
## 1656.4026 1679.2414 -821.2013       186 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups   Name        Std.Dev.
##  Batch    (Intercept)  4.633  
##  Residual             16.482  
## 
## Number of obs: 193 / Conditional model: Batch, 50
## 
## Dispersion estimate for gaussian family (sigma^2):  272 
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)     DeathAge    IsolateC1   IsolateC18   IsolateC24  
##    -20.3723       1.3133      -4.2733      -0.6133      -6.5511  
## 
## $<NA>
## NULL
## 
## attr(,"rank")
## function (x) 
## do.call("rank", list(x))
## <environment: 0x137323810>
## attr(,"call")
## AICc(x)
## attr(,"class")
## [1] "function"     "rankFunction"
## attr(,"beta")
## [1] "none"
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.78 0.656 0.584 0.568 0.468 0.564 0.596 0.448 0.512 0.388 0.464 0.368 0.504 0.452 0.672 0.536 0.708 0.592 0.356 0.488 ...
tab_model(best_model, show.intercept = F) # No effect
  Nb_offspring
Predictors Estimates CI p
DeathAge 1.31 1.16 – 1.45 <0.001
Random Effects
σ2 279.59
τ00 Batch 21.15
ICC 0.07
N Batch 50
Observations 193
Marginal R2 / Conditional R2 0.622 / 0.649 
summary(best_model)
##  Family: gaussian  ( identity )
## Formula:          Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
## 
##      AIC      BIC   logLik deviance df.resid 
##   1655.5   1668.5   -823.7   1647.5      189 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  Batch    (Intercept)  21.15    4.599  
##  Residual             279.59   16.721  
## Number of obs: 193, groups:  Batch, 50
## 
## Dispersion estimate for gaussian family (sigma^2):  280 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -22.63107    4.34365   -5.21 1.89e-07 ***
## DeathAge      1.30586    0.07379   17.70  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculation of daily reproduction rate
repro_nf$daily_repro_rate = repro_nf$Nb_offspring / repro_nf$DeathAge
mean(repro_nf$daily_repro_rate, na.rm = TRUE)
## [1] 0.8219201
(sem = sd(repro_nf$daily_repro_rate, na.rm = TRUE) / sqrt(length(repro_nf$daily_repro_rate)))
## [1] 0.02774078
Visualisation

Only the effect of time is important here; each day, a Daphnia produces on average 0.30 more offspring.

4.2.3 Uninfected individuals (including unexposed husbandry control)

dataC$Isolate = as.character(dataC$Isolate)
dataC$Isolate[is.na(dataC$Isolate)] = "unexposed"
dataC$Isolate = as.factor(dataC$Isolate)

dataC$Temperature = as.character(dataC$Temperature)
dataC$Temperature[dataC$Temperature == "CTRL"] = 20
dataC$Temperature = factor(dataC$Temperature, ordered = TRUE)

repro_nfC = subset(dataC, dataC$Clutch>0 & InfectedStatus!="inf")

repro_nfC$Temperature = factor(repro_nfC$Temperature, labels = c("20", "30","40"), ordered = T)

FC1 = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature+ DeathAge, family = "poisson", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)

simulateResiduals(FC1, plot=T) # I don't find better and I have no reason to exclude the outlier.

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.724 0.608 0.58 0.488 0.6 0.48 0.656 0.54 0.672 0.492 0.484 0.408 0.4 0.548 0.464 0.352 0.8 0.624 0.824 0.668 ...
FC1R = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge  + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nfC, family = "gaussian", REML=F)

AICc(FC1, FC1R)# Batch effect NOT important and do not improve residuals anyway.
##      df     AICc
## FC1   9 1945.748
## FC1R 10 1946.820
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge, data = repro_nfC, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1) 

model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge, 
##     data = repro_nfC, family = "gaussian", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df    logLik
## 2      -22.74          +    1.305                                 3  -966.679
## 4      -24.62          +    1.311        +                        7  -963.982
## 6      -22.39          +    1.295                 +               5  -966.198
## 8      -24.57          +    1.305        +        +               9  -963.461
## 16     -24.28          +    1.314        +        +            + 15  -962.257
## 5       47.93          +                          +               4 -1078.562
## 1       48.49          +                                          2 -1083.107
## 7       50.74          +                 +        +               8 -1077.345
## 3       52.70          +                 +                        6 -1079.491
## 15      46.80          +                 +        +            + 14 -1076.437
##      AICc  delta weight
## 2  1939.5   0.00  0.682
## 4  1942.5   3.01  0.151
## 6  1942.7   3.20  0.137
## 8  1945.7   6.28  0.029
## 16 1956.8  17.31  0.000
## 5  2165.3 225.84  0.000
## 1  2170.3 230.80  0.000
## 7  2171.3 231.88  0.000
## 3  2171.4 231.90  0.000
## 15 2182.8 243.38  0.000
## Models ranked by AICc(x)
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.776 0.628 0.604 0.564 0.528 0.548 0.648 0.632 0.548 0.436 0.528 0.432 0.496 0.436 0.432 0.452 0.724 0.516 0.784 0.664 ...
selected_models = get.models(model_selection, 1:nrow(model_selection))
tab_model(selected_models, show.intercept = F) # No effect
  Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
DeathAge 1.30 1.18 – 1.43 <0.001 1.31 1.18 – 1.44 <0.001 1.29 1.17 – 1.42 <0.001 1.31 1.18 – 1.43 <0.001 1.31 1.18 – 1.44 <0.001
Isolate [C14] 4.03 -2.97 – 11.04 0.259 3.83 -3.22 – 10.88 0.287 3.21 -4.95 – 11.37 0.441 1.19 -10.42 – 12.80 0.841 0.82 -10.79 – 12.43 0.890 5.39 -8.07 – 18.85 0.433
Isolate [C18] 3.72 -3.09 – 10.52 0.284 3.82 -3.00 – 10.63 0.272 3.17 -4.64 – 10.98 0.426 -3.36 -14.52 – 7.80 0.555 -4.35 -15.56 – 6.87 0.447 0.98 -11.90 – 13.86 0.882
Isolate [C24] -2.42 -9.45 – 4.62 0.501 -2.51 -9.64 – 4.63 0.491 -2.96 -11.08 – 5.16 0.475 -5.37 -17.12 – 6.37 0.370 -6.47 -18.12 – 5.19 0.277 -1.23 -14.62 – 12.17 0.858
Isolate [CTRL] 1.21 -7.28 – 9.69 0.780 2.25 -7.51 – 12.00 0.652 -0.81 -18.17 – 16.54 0.927 -7.95 -23.93 – 8.03 0.330 -15.55 -29.36 – -1.74 0.027 7.65 -20.95 – 36.26 0.600
Temperature [linear] 1.63 -2.10 – 5.37 0.391 1.71 -2.54 – 5.96 0.431 -1.60 -14.45 – 11.26 0.808 9.22 3.23 – 15.20 0.003 7.35 0.41 – 14.29 0.038 19.74 -1.19 – 40.67 0.065
Temperature [quadratic] 0.81 -3.19 – 4.80 0.692 0.98 -3.23 – 5.18 0.649 1.01 -8.88 – 10.90 0.841 -2.19 -8.71 – 4.33 0.510 -1.09 -8.01 – 5.83 0.757 -8.20 -24.46 – 8.05 0.323
Isolate [C14] ×
Temperature [linear] 		3.67 -11.49 – 18.83 0.635 -14.13 -38.99 – 10.72 0.265
Isolate [C18] ×
Temperature [linear] 		3.79 -11.12 – 18.70 0.618 -14.78 -39.20 – 9.63 0.235
Isolate [C24] ×
Temperature [linear] 		3.84 -10.96 – 18.64 0.611 -12.61 -36.88 – 11.67 0.309
Isolate [C14] ×
Temperature [quadratic] 		-2.63 -15.81 – 10.55 0.696 8.53 -13.14 – 30.20 0.440
Isolate [C18] ×
Temperature [quadratic] 		4.00 -8.12 – 16.11 0.518 9.27 -10.71 – 29.24 0.363
Isolate [C24] ×
Temperature [quadratic] 		-2.94 -16.34 – 10.46 0.667 5.90 -16.16 – 27.96 0.600
Observations 228 228 228 228 228 228 228 228 228 228
R2 / R2 adjusted 0.640 / 0.637 0.648 / 0.639 0.641 / 0.635 0.650 / 0.637 0.654 / 0.631 0.039 / 0.026 0.000 / -0.004 0.049 / 0.019 0.031 / 0.009 0.057 / -0.000 
summary(best_model)
##  Family: gaussian  ( identity )
## Formula:          Nb_offspring ~ DeathAge + 1
## Data: repro_nfC
## 
##      AIC      BIC   logLik deviance df.resid 
##   1939.4   1949.6   -966.7   1933.4      225 
## 
## 
## Dispersion estimate for gaussian family (sigma^2):  282 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -22.74119    3.70948  -6.131 8.76e-10 ***
## DeathAge      1.30473    0.06482  20.127  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Again, the husbandry control do not change the conclusions.

Visualisation

4.2.4 Descriptive stat per group

Here we just described the average number of offspring produced per day in total in each group

# Function to compute mean and 95% confidence interval
mean_ci = function(x) {
  x = x[!is.na(x)]
  m = mean(x)
  se = sd(x) / sqrt(length(x))
  c(mean = m,
    CI_low = m - 1.96 * se,
    CI_high = m + 1.96 * se)
}

# Define groups
g1 = mean_ci(repro_nfC$repro_rate[repro_nfC$Isolate == "CTRL"])
g2 = mean_ci(repro_nf$repro_rate[repro_nfC$Temperature == 20])
g3 = mean_ci(repro_nf$repro_rate)
g4 = mean_ci(repro_inf$repro_rate)

rbind(
  noninfected_unexposed     = g1,
  noninfected_exposed_20    = g2,
  noninfected_exposed_all   = g3,
  infected_all              = g4
)
##                              mean    CI_low   CI_high
## noninfected_unexposed   0.7203259 0.5747725 0.8658793
## noninfected_exposed_20  0.7463513 0.6490932 0.8436094
## noninfected_exposed_all 0.8219201 0.7675482 0.8762921
## infected_all            0.1699311 0.1302886 0.2095736

5. Parasite life-cycle

This section explores how temperature and the genotype of the parasite contribute to the abundance of each stage of the life-cycle after infection, using generalized linear mixed models.

5.1 Cauliflowers

Cauliflowers corresponds to the first stage of the parasite life-cycle and here we are interested to understand how their presence/absence is affected by the temperature and/or the genotype of the spores. We also investigated the potential impact of the variations in time post-exposure, as all screened individuals died at different time.

data_c1 = subset(data, data$cauliflowers != "not-well-counted" & data$cauliflowers != "counted-after-freezer" & data$stade_max != "2nd_Cauliflower" & Temperature != "NA")
data_c1$Temperature = factor(data_c1$Temperature, ordered = TRUE)

Caul = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate+DeathAge, family = "binomial", REML= T)
CaulR = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate+DeathAge + (1|Batch), family = "binomial", REML= T)

AICc(Caul, CaulR)
##       df     AICc
## Caul   7 119.6534
## CaulR  8 121.7682
Caul1 = glmmTMB(data=data_c1, cauliflowers~Temperature*Isolate*DeathAge, family = "binomial", REML= F)
Caul2 = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate*DeathAge, family = "binomial", REML= F)
Caul3 = glmmTMB(data=data_c1, cauliflowers~Temperature*Isolate+DeathAge, family = "binomial", REML= F) 
Caul4 = glmmTMB(data=data_c1, cauliflowers~Isolate+DeathAge*Temperature, family = "binomial", REML= F)
Caul5 = glmmTMB(data=data_c1, cauliflowers~Isolate*DeathAge, family = "binomial", REML= F)
Caul6 = glmmTMB(data=data_c1, cauliflowers~Temperature*Isolate, family = "binomial", REML= F)
Caul7 = glmmTMB(data=data_c1, cauliflowers~DeathAge*Temperature, family = "binomial", REML= F)
Caul8 = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate+DeathAge, family = "binomial", REML= F)
Caul9 = glmmTMB(data=data_c1, cauliflowers~Isolate+DeathAge, family = "binomial", REML= F)
Caul10 = glmmTMB(data=data_c1, cauliflowers~Temperature+Isolate, family = "binomial", REML= F)
Caul11 = glmmTMB(data=data_c1, cauliflowers~DeathAge+Temperature, family = "binomial", REML= F)
Caul12 = glmmTMB(data=data_c1, cauliflowers~DeathAge, family = "binomial", REML= F)
Caul13 = glmmTMB(data=data_c1, cauliflowers~Isolate, family = "binomial", REML= F)
Caul14 = glmmTMB(data=data_c1, cauliflowers~Temperature, family = "binomial", REML= F)
Caul0 = glmmTMB(data=data_c1, cauliflowers~1, family = "binomial", REML= F)

AICc(Caul1, Caul2, Caul3, Caul4, Caul5, Caul6, Caul7, Caul8, Caul9, Caul10, Caul11, Caul12, Caul13, Caul14, Caul0)
##        df     AICc
## Caul1  24       NA
## Caul2  10 121.1953
## Caul3  13 123.4742
## Caul4   9 121.1437
## Caul5   8 121.2525
## Caul6  12 151.2119
## Caul7   6 115.2956
## Caul8   7 116.9950
## Caul9   5 117.3381
## Caul10  6 148.7161
## Caul11  4 111.2346
## Caul12  2 111.9709
## Caul13  4 157.5159
## Caul14  3 144.1507
## Caul0   1 152.8702
tab_model(Caul11, Caul12, Caul7, Caul8, Caul5, Caul2, Caul4, show.intercept = F)
  cauliflowers cauliflowers cauliflowers cauliflowers cauliflowers cauliflowers cauliflowers
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
DeathAge 0.91 0.87 – 0.95 <0.001 0.90 0.87 – 0.94 <0.001 0.91 0.87 – 0.95 <0.001 0.91 0.87 – 0.95 <0.001 0.85 0.73 – 1.00 0.049 0.86 0.73 – 1.01 0.064 0.91 0.87 – 0.95 <0.001
Temperature [linear] 0.34 0.11 – 1.05 0.060 0.24 0.02 – 3.09 0.271 0.34 0.11 – 1.08 0.066 0.34 0.10 – 1.08 0.068 0.23 0.02 – 3.33 0.284
Temperature [quadratic] 0.97 0.36 – 2.57 0.944 0.84 0.08 – 9.24 0.886 0.96 0.36 – 2.57 0.940 0.94 0.35 – 2.51 0.895 0.81 0.07 – 9.22 0.866
DeathAge × Temperature
[linear] 		1.01 0.93 – 1.10 0.753 1.01 0.93 – 1.10 0.755
DeathAge × Temperature
[quadratic] 		1.00 0.92 – 1.09 0.911 1.01 0.92 – 1.10 0.890
Isolate [C1] 0.91 0.21 – 3.99 0.897 0.16 0.00 – 12.00 0.409 0.11 0.00 – 9.38 0.331 0.89 0.20 – 3.94 0.878
Isolate [C18] 0.85 0.19 – 3.83 0.831 0.06 0.00 – 4.98 0.214 0.10 0.00 – 10.12 0.331 0.87 0.19 – 3.99 0.857
Isolate [C24] 1.33 0.34 – 5.16 0.683 0.91 0.01 – 134.04 0.969 1.02 0.01 – 159.66 0.994 1.34 0.34 – 5.23 0.676
Isolate [C1] × DeathAge 1.08 0.91 – 1.29 0.389 1.09 0.91 – 1.31 0.340
Isolate [C18] × DeathAge 1.10 0.92 – 1.31 0.285 1.09 0.91 – 1.30 0.363
Isolate [C24] × DeathAge 1.02 0.82 – 1.25 0.886 1.01 0.82 – 1.25 0.899
Observations 291 291 291 291 291 291 291
R2 Tjur 0.194 0.154 0.197 0.203 0.183 0.220 0.205 
# No differences

simulateResiduals(Caul11, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3388901 0.3332647 0.415397 0.04419869 0.3642384 0.5669417 0.7992067 0.6843859 0.5363047 0.5255567 0.4620051 0.9983571 0.2413181 0.1928017 0.1345064 0.2465995 0.7431666 0.4926976 0.5356813 0.0344749 ...
# best model 
Caul11 = glmmTMB(data=data_c1, cauliflowers~DeathAge+Temperature, family = "binomial", REML=T)
summary(Caul11)
##  Family: binomial  ( logit )
## Formula:          cauliflowers ~ DeathAge + Temperature
## Data: data_c1
## 
##      AIC      BIC   logLik deviance df.resid 
##    116.8    131.4    -54.4    108.8      287 
## 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.66549    0.67120   0.991   0.3214    
## DeathAge      -0.09664    0.02279  -4.241 2.23e-05 ***
## Temperature.L -1.08254    0.57628  -1.878   0.0603 .  
## Temperature.Q -0.03522    0.49967  -0.070   0.9438    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(Caul11, show.intercept = F)
  cauliflowers
Predictors Odds Ratios CI p
DeathAge 0.91 0.87 – 0.95 <0.001
Temperature [linear] 0.34 0.11 – 1.05 0.060
Temperature [quadratic] 0.97 0.36 – 2.57 0.944
Observations 291
R2 Tjur 0.194


We observed that the time post-exposure is reducing of 10% the probability to observe first stage cauliflowers, temperature might have a role that is not significant with the current effective of our experimentation.

Visualisation

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_line()`).

5.2 Grape seed stage

The grape seed stage corresponds to the second visible stage of the parasite life cycle and here we are interested in understanding how its presence/absence is affected by temperature and the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.

data_crushed_inf_c = subset(infected, infected$grapes!="counted-after-freezer" & infected$grapes!="counted-after-freezer")

Grap = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate+DeathAge, family = "binomial", REML=F)
GrapR = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate+DeathAge+ (1|Batch), family = "binomial", REML=F)

AICc(Grap,GrapR) # No need for random effect
##       df     AICc
## Grap   7 209.1698
## GrapR  8 210.7109
Grap1 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*Isolate*DeathAge, family = "binomial")
Grap2 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate*DeathAge, family = "binomial")
Grap3 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*Isolate+DeathAge, family = "binomial")
Grap4 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate+DeathAge*Temperature, family = "binomial")

Grap5 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate*DeathAge, family = "binomial")
Grap6 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*Isolate, family = "binomial")
Grap7 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature*DeathAge, family = "binomial")

Grap8 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate+DeathAge, family = "binomial")

Grap9 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate+DeathAge, family = "binomial")
Grap10 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+Isolate, family = "binomial")
Grap11 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature+DeathAge, family = "binomial")

Grap12 = glmmTMB(data=data_crushed_inf_c, grapes~DeathAge, family = "binomial")
Grap13 = glmmTMB(data=data_crushed_inf_c, grapes~Isolate, family = "binomial")
Grap14 = glmmTMB(data=data_crushed_inf_c, grapes~Temperature, family = "binomial")

Grap0 = glmmTMB(data=data_crushed_inf_c, grapes~1, family = "binomial")

AICc(Grap1, Grap2, Grap3, Grap4, Grap5, Grap6, Grap7, Grap8, Grap9, Grap10, Grap11, Grap12, Grap13, Grap14, Grap0)
##        df     AICc
## Grap1  24       NA
## Grap2  10 204.7698
## Grap3  13 219.5599
## Grap4   9 212.0557
## Grap5   8 207.6257
## Grap6  12 217.6432
## Grap7   6 206.7304
## Grap8   7 209.1698
## Grap9   5 213.9978
## Grap10  6 207.8652
## Grap11  4 204.1768
## Grap12  2 209.5133
## Grap13  4 213.0738
## Grap14  3 202.8924
## Grap0   1 208.6435
tab_model(Grap14, Grap11, Grap10, Grap8, Grap0, show.intercept = F)
  grapes grapes grapes grapes grapes
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
Temperature [linear] 2.69 1.38 – 5.25 0.004 2.66 1.36 – 5.20 0.004 2.70 1.37 – 5.32 0.004 2.67 1.35 – 5.27 0.005
Temperature [quadratic] 1.44 0.79 – 2.62 0.237 1.41 0.77 – 2.57 0.268 1.40 0.76 – 2.56 0.278 1.37 0.74 – 2.51 0.315
DeathAge 1.01 0.99 – 1.04 0.366 1.01 0.99 – 1.04 0.348
Isolate [C1] 1.76 0.63 – 4.89 0.281 1.84 0.65 – 5.18 0.249
Isolate [C18] 1.57 0.59 – 4.17 0.363 1.44 0.53 – 3.90 0.478
Isolate [C24] 1.25 0.47 – 3.36 0.652 1.22 0.45 – 3.29 0.693
Observations 158 158 158 158 158
R2 Tjur 0.058 0.064 0.068 0.075 0.000 
# Temperature increase the chances to observe grapes

simulateResiduals(Grap14, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.7676696 0.05811451 0.2656332 0.3575351 0.9665073 0.3685181 0.7259261 0.04275554 0.06533833 0.1598172 0.02149855 0.8880048 0.4096685 0.3173512 0.08902606 0.6813057 0.1156503 0.1831655 0.3961829 0.6764339 ...
summary(Grap14)
##  Family: binomial  ( logit )
## Formula:          grapes ~ Temperature
## Data: data_crushed_inf_c
## 
##      AIC      BIC   logLik deviance df.resid 
##    202.7    211.9    -98.4    196.7      155 
## 
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.7337     0.1873   3.918 8.94e-05 ***
## Temperature.L   0.9886     0.3416   2.894   0.0038 ** 
## Temperature.Q   0.3618     0.3062   1.182   0.2373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(Grap14, show.intercept = F)
  grapes
Predictors Odds Ratios CI p
Temperature [linear] 2.69 1.38 – 5.25 0.004
Temperature [quadratic] 1.44 0.79 – 2.62 0.237
Observations 158
R2 Tjur 0.058


We observed that the temperature increases the odds of observing grape stage seeds by 168%.

Visualisation

5.3 Immature spores

After grape seed formation, spores begin their final maturation stage and develop their final shape. However, this is gradual and when an individual died we could observe partially matured spores, which we called immature spores. Here we are interested in understanding how their presence/absence is affected by temperature and/or the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.

data_crushed_inf_c = subset(infected, infected$imm_load!="counted-after-freezer")
data_crushed_inf_c$imm_presence = ifelse(data_crushed_inf_c$imm >0, 1,0)

ImmP = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate+DeathAge, family = "binomial")
ImmPR = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate+DeathAge+ (1|Batch), family = "binomial")

AICc(ImmP,ImmPR) # No need for random effect
##       df     AICc
## ImmP   7 109.2572
## ImmPR  8 108.0838
ImmP1 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*Isolate*DeathAge, family = "binomial")
ImmP2 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate*DeathAge, family = "binomial")
ImmP3 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*Isolate+DeathAge, family = "binomial")
ImmP4 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate+DeathAge*Temperature, family = "binomial")

ImmP5 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate*DeathAge, family = "binomial")
ImmP6 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*Isolate, family = "binomial")
ImmP7 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature*DeathAge, family = "binomial")

ImmP8 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate+DeathAge, family = "binomial")

ImmP9 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate+DeathAge, family = "binomial")
ImmP10 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+Isolate, family = "binomial")
ImmP11 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature+DeathAge, family = "binomial")

ImmP12 = glmmTMB(data=data_crushed_inf_c, imm_presence~DeathAge, family = "binomial")
ImmP13 = glmmTMB(data=data_crushed_inf_c, imm_presence~Isolate, family = "binomial")
ImmP14 = glmmTMB(data=data_crushed_inf_c, imm_presence~Temperature, family = "binomial")

ImmP0 = glmmTMB(data=data_crushed_inf_c, imm_presence~1, family = "binomial")

AICc(ImmP1, ImmP2, ImmP3, ImmP4, ImmP5, ImmP6, ImmP7, ImmP8, ImmP9, ImmP10, ImmP11, ImmP12, ImmP13, ImmP14, ImmP0)
##        df     AICc
## ImmP1  24       NA
## ImmP2  10 111.6544
## ImmP3  13 112.0329
## ImmP4   9 110.5120
## ImmP5   8 110.1232
## ImmP6  12 146.8456
## ImmP7   6 104.4657
## ImmP8   7 109.2572
## ImmP9   5 108.3799
## ImmP10  6 148.4499
## ImmP11  4 103.2475
## ImmP12  2 102.5769
## ImmP13  4 147.5806
## ImmP14  3 143.6712
## ImmP0   1 143.0355
tab_model(ImmP12, ImmP11, ImmP7, show.intercept = F)
  imm_presence imm_presence imm_presence
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
DeathAge 1.13 1.08 – 1.19 <0.001 1.14 1.08 – 1.19 <0.001 1.13 1.07 – 1.20 <0.001
Temperature [linear] 1.20 0.48 – 3.00 0.690 1.48 0.03 – 80.40 0.849
Temperature [quadratic] 0.43 0.15 – 1.20 0.106 0.04 0.00 – 0.99 0.050
Temperature [linear] ×
DeathAge 		0.99 0.89 – 1.10 0.882
Temperature [quadratic] ×
DeathAge 		1.08 0.99 – 1.17 0.101
Observations 167 167 167
R2 Tjur 0.311 0.344 0.361 
# Temperature increase the chances to observe imm_presence

simulateResiduals(ImmP12, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.08681092 0.4781304 0.5053635 0.3460266 0.9493186 0.5617266 0.07295504 0.2366983 0.1756649 0.6546839 0.7908228 0.8876897 0.893821 0.7406902 0.239505 0.7249811 0.1957589 0.5392972 0.5591017 0.8727609 ...
summary(ImmP12)
##  Family: binomial  ( logit )
## Formula:          imm_presence ~ DeathAge
## Data: data_crushed_inf_c
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.5    108.7    -49.3     98.5      165 
## 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.12919    0.84637  -3.697 0.000218 ***
## DeathAge     0.12522    0.02326   5.384 7.27e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(ImmP12, show.intercept = F)
  imm_presence
Predictors Odds Ratios CI p
DeathAge 1.13 1.08 – 1.19 <0.001
Observations 167
R2 Tjur 0.311


We observed that only the time between exposure and death is increase of 13% the probability to observe immmature presence after host death, which is why we removed this abalysis from the main manuscript.

Visualisation

5.4 Mature spores

Spores reach their final maturation stage and can now infect new hosts. We are interested in understanding how individual mature parasite load is affected by temperature and/or the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.

data_mat = subset(infected, infected$mature_load>0)
data_mat$Temperature = factor(data_mat$Temperature, ordered = F)

mat1 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate*DeathAge, family = "gaussian", REML=T)
mat1R = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate*DeathAge+(1|Batch), family = "gaussian", REML=T)

AICc(mat1, mat1R) # No need of a batch effect
##       df     AICc
## mat1  25 242.8509
## mat1R 26 245.6588
mat2 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+Isolate*DeathAge, family = "gaussian")
mat3 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate+DeathAge, family = "gaussian")
mat4 = glmmTMB(data=data_mat, mature_load/10000000~Isolate+DeathAge*Temperature, family = "gaussian")

mat5 = glmmTMB(data=data_mat, mature_load/10000000~Isolate*DeathAge, family = "gaussian")
mat6 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*Isolate, family = "gaussian")
mat7 = glmmTMB(data=data_mat, mature_load/10000000~Temperature*DeathAge, family = "gaussian")

mat8 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+Isolate+DeathAge, family = "gaussian")

mat9 = glmmTMB(data=data_mat, mature_load/10000000~Isolate+DeathAge, family = "gaussian")
mat10 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+Isolate, family = "gaussian")
mat11 = glmmTMB(data=data_mat, mature_load/10000000~Temperature+DeathAge, family = "gaussian")

mat12 = glmmTMB(data=data_mat, mature_load/10000000~DeathAge, family = "gaussian")
mat13 = glmmTMB(data=data_mat, mature_load/10000000~Isolate, family = "gaussian")
mat14 = glmmTMB(data=data_mat, mature_load/10000000~Temperature, family = "gaussian")

mat0 = glmmTMB(data=data_mat, mature_load/10000000~1, family = "gaussian")

AICc(mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10, mat11, mat12, mat13, mat14, mat0)
##       df      AICc
## mat1  25 242.85092
## mat2  11  97.42244
## mat3  14 104.23161
## mat4  10  99.48683
## mat5   9  93.04474
## mat6  13 150.21483
## mat7   7 107.51440
## mat8   8  96.60519
## mat9   6  92.28979
## mat10  7 139.42053
## mat11  5 104.35369
## mat12  3 100.24447
## mat13  5 136.94938
## mat14  4 154.88560
## mat0   2 152.62650
tab_model(mat9, mat5, show.intercept = F)
  mature_load/1e+07 mature_load/1e+07
Predictors Estimates CI p Estimates CI p
Isolate [C1] 0.09 -0.07 – 0.24 0.286 -0.70 -1.46 – 0.06 0.071
Isolate [C18] 0.28 0.13 – 0.43 <0.001 -0.44 -1.10 – 0.22 0.189
Isolate [C24] 0.09 -0.07 – 0.24 0.286 -0.35 -1.03 – 0.33 0.308
DeathAge 0.02 0.01 – 0.02 <0.001 0.01 -0.00 – 0.02 0.122
Isolate [C1] × DeathAge 0.02 0.00 – 0.03 0.043
Isolate [C18] × DeathAge 0.01 0.00 – 0.03 0.028
Isolate [C24] × DeathAge 0.01 -0.00 – 0.02 0.184
Observations 145 145
R2 / R2 adjusted 0.378 / 0.356 0.403 / 0.368 
simulateResiduals(mat5, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.468 0.368 0.34 0.828 0.316 0.808 0.128 0.4 0.48 0.48 0.124 0.096 0.224 0.332 0.36 0.516 0.54 0.3 0.34 0.48 ...
tab_model(mat5, show.intercept = F)
  mature_load/1e+07
Predictors Estimates CI p
Isolate [C1] -0.70 -1.46 – 0.06 0.071
Isolate [C18] -0.44 -1.10 – 0.22 0.189
Isolate [C24] -0.35 -1.03 – 0.33 0.308
DeathAge 0.01 -0.00 – 0.02 0.122
Isolate [C1] × DeathAge 0.02 0.00 – 0.03 0.043
Isolate [C18] × DeathAge 0.01 0.00 – 0.03 0.028
Isolate [C24] × DeathAge 0.01 -0.00 – 0.02 0.184
Observations 145
R2 / R2 adjusted 0.403 / 0.368


C1 and C18 have significantly increased growth speed over time compared to C14 (référence) and C24.

Visualisation

5.5 2nd cauliflowers

Cauliflowers can also be observed if parasites restart a second cycle inside the same host. Here we are interested in understanding how their presence/absence is affected by temperature and the genotype of the spores. We also investigated the potential impact of variations in time post-exposure, as all screened individuals died at different times.

data_c2 = subset(infected, infected$cauliflowers!="counted-after-freezer" & infected$cauliflowers!="counted-after-freezer"& infected$stade_max!="Cauliflower")

data_c2$Temperature = factor(data_c2$Temperature, ordered = T)

Caul2 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate+DeathAge, family = "binomial")
Caul2R = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate+DeathAge + (1|Batch), family = "binomial", REML= T)

AICc(Caul2, Caul2R)
##        df     AICc
## Caul2  13 207.2512
## Caul2R 14 207.1492
Caul2.1 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate*DeathAge, family = "binomial")

Caul2.2 = glmmTMB(data=data_c2, cauliflowers~Temperature+Isolate*DeathAge, family = "binomial")
Caul2.3 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate+DeathAge, family = "binomial")
Caul2.4 = glmmTMB(data=data_c2, cauliflowers~Isolate+DeathAge*Temperature, family = "binomial")

Caul2.5 = glmmTMB(data=data_c2, cauliflowers~Isolate*DeathAge, family = "binomial")
Caul2.6 = glmmTMB(data=data_c2, cauliflowers~Temperature*Isolate, family = "binomial")
Caul2.7 = glmmTMB(data=data_c2, cauliflowers~DeathAge*Temperature, family = "binomial")

Caul2.8 = glmmTMB(data=data_c2, cauliflowers~Temperature+Isolate+DeathAge, family = "binomial")

Caul2.9 = glmmTMB(data=data_c2, cauliflowers~Isolate+DeathAge, family = "binomial")
Caul2.10 = glmmTMB(data=data_c2, cauliflowers~Temperature+Isolate, family = "binomial")
Caul2.11 = glmmTMB(data=data_c2, cauliflowers~DeathAge+Temperature, family = "binomial")

Caul2.12 = glmmTMB(data=data_c2, cauliflowers~DeathAge, family = "binomial")
Caul2.13 = glmmTMB(data=data_c2, cauliflowers~Isolate, family = "binomial")
Caul2.14 = glmmTMB(data=data_c2, cauliflowers~Temperature, family = "binomial")

Caul2.0 = glmmTMB(data=data_c2, cauliflowers~1, family = "binomial")

AICc(Caul2.1, Caul2.2, Caul2.3, Caul2.4, Caul2.5, Caul2.6, Caul2.7, Caul2.8, Caul2.9, Caul2.10, Caul2.11, Caul2.12, Caul2.13, Caul2.14, Caul2.0)
##          df     AICc
## Caul2.1  24 221.5962
## Caul2.2  10 206.3934
## Caul2.3  13 207.2512
## Caul2.4   9 206.1231
## Caul2.5   8 205.8774
## Caul2.6  12 205.5078
## Caul2.7   6 205.5147
## Caul2.8   7 201.9549
## Caul2.9   5 201.6345
## Caul2.10  6 201.1274
## Caul2.11  4 201.8077
## Caul2.12  2 201.3683
## Caul2.13  4 200.2414
## Caul2.14  3 199.9651
## Caul2.0   1 199.3769
tab_model(Caul2.0, Caul2.14, Caul2.13, Caul2.12, Caul2.10, show.intercept = F)
  cauliflowers cauliflowers cauliflowers cauliflowers cauliflowers
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
Temperature [linear] 1.06 0.60 – 1.89 0.837 1.04 0.58 – 1.89 0.891
Temperature [quadratic] 0.58 0.32 – 1.05 0.070 0.57 0.31 – 1.05 0.073
Isolate [C1] 0.91 0.32 – 2.61 0.858 0.98 0.33 – 2.85 0.965
Isolate [C18] 0.38 0.14 – 1.02 0.054 0.39 0.14 – 1.07 0.067
Isolate [C24] 0.69 0.25 – 1.94 0.488 0.72 0.25 – 2.05 0.535
DeathAge 1.00 0.98 – 1.03 0.799
Observations 145 145 145 145 145
R2 Tjur 0.000 0.024 0.037 0.000 0.059 
# No differences

simulateResiduals(Caul2.13, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.2323121 0.4716365 0.7126411 0.6861518 0.2546818 0.3056335 0.645681 0.6405459 0.9225231 0.8999696 0.9554663 0.1963001 0.20196 0.6420615 0.6905008 0.00794358 0.1867782 0.665696 0.4064166 0.2371338 ...
tab_model(Caul2.10, show.intercept = F)
  cauliflowers
Predictors Odds Ratios CI p
Temperature [linear] 1.04 0.58 – 1.89 0.891
Temperature [quadratic] 0.57 0.31 – 1.05 0.073
Isolate [C1] 0.98 0.33 – 2.85 0.965
Isolate [C18] 0.39 0.14 – 1.07 0.067
Isolate [C24] 0.72 0.25 – 2.05 0.535
Observations 145
R2 Tjur 0.059


We observed that the null model is among the best ones and that the measured effect are relatively low, so we preffered the model taking into account the effect of the genotype on the second cauliflower formation, which seem more biologically relevant, however it should be interptreted with caution.

Visualisation


The isolate C18 is a little bit less likely to start a second cycle inside it’s host thant the other lineages (66% less in average).

5.6 Maximal stage reached

The purpose of this section is to analyze the impact of temperature and parasite genotype on the maximum stage reached by the parasite at host death. We also included the impact of host age at death as individuals died at different times post-exposure.

infected = subset(infected,stade_max!="None")
infected$stade_max = factor(infected$stade_max, ordered = T)
infected$Temperature = factor(infected$Temperature, ordered = T)

# On crée une nouvelle variable pour ne pas écraser l'originale
infected$stade_simplifie <- fct_collapse(infected$stade_max,
  "Caul_Grapes" = c("Cauliflower", "Grapes"),
  "Imm_Mature"  = c("Immature", "Mature"),
  "2nd_Caul"    = "2nd_Cauliflower"
)

# Crucial : On s'assure que R comprend que c'est ordonné
# (Vérifie bien que cet ordre est le bon ordre biologique !)
infected$stade_simplifie <- factor(infected$stade_simplifie, 
                                  levels = c("Caul_Grapes", "Imm_Mature", "2nd_Caul"), 
                                  ordered = TRUE)

#modele13 = polr(stade_max ~ Temperature * Isolate * DeathAge, data = infected, Hess = TRUE)

modele14 = polr(stade_max ~ Temperature + Isolate + DeathAge, data = infected, Hess = TRUE)
modele14R = clmm(stade_max ~ Temperature + Isolate + DeathAge + (1 | Batch), data = infected)
AICc(modele14R, modele14) # No need for random factor
##           df     AICc
## modele14R 11 293.3242
## modele14  10 292.0298
#modele1 = polr(stade_max ~ Temperature + Isolate * DeathAge, data = infected, Hess = TRUE)
modele2 = polr(stade_max ~ Temperature * Isolate + DeathAge, data = infected, Hess = TRUE)
modele3 = polr(stade_max ~ Isolate +  DeathAge * Temperature, data = infected, Hess = TRUE)

#modele4 = polr(stade_max ~ Isolate * DeathAge, data = infected, Hess = TRUE)
modele5 = polr(stade_max ~ Temperature * Isolate, data = infected, Hess = TRUE)
modele6 = polr(stade_max ~ DeathAge * Temperature, data = infected, Hess = TRUE)

modele7 = polr(stade_max ~ Isolate + DeathAge, data = infected, Hess = TRUE)
modele8 = polr(stade_max ~ Temperature + Isolate, data = infected, Hess = TRUE)
modele9 = polr(stade_max ~ DeathAge + Temperature, data = infected, Hess = TRUE)

modele10 = polr(stade_max ~ Temperature, data = infected, Hess = TRUE)
modele11 = polr(stade_max ~ Isolate, data = infected, Hess = TRUE)
modele12 = polr(stade_max ~ DeathAge, data = infected, Hess = TRUE)

modele0 = polr(stade_max ~ 1, data = infected, Hess = TRUE)

AICc(modele2,modele3,modele5,modele6,modele7,modele8,modele9,modele10,modele11,modele12,modele14,modele0)
##          df     AICc
## modele2  16 301.0434
## modele3  12 293.9404
## modele5  15 347.5868
## modele6   9 297.1830
## modele7   8 299.5568
## modele8   9 341.0864
## modele9   7 294.9343
## modele10  6 335.5819
## modele11  7 342.4234
## modele12  5 301.1638
## modele14 10 292.0298
## modele0   4 337.1707
tab_model(modele14, modele3, modele9)
  stade_max stade_max stade_max
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
Cauliflower|Grapes 3.24 0.85 – 12.35 0.084 2.01 0.44 – 9.15 0.364 2.60 0.81 – 8.37 0.109
Grapes|Immature 3.60 0.94 – 13.75 0.061 2.24 0.49 – 10.19 0.294 2.89 0.89 – 9.31 0.076
Immature|Mature 4.40 1.14 – 16.94 0.032 2.75 0.60 – 12.52 0.190 3.52 1.08 – 11.43 0.037
Mature|2nd Cauliflower 95.60 19.81 – 461.41 <0.001 61.98 11.36 – 338.24 <0.001 63.92 16.10 – 253.77 <0.001
Temperature [linear] 1.53 0.87 – 2.73 0.145 9.54 0.99 – 104.48 0.057 1.47 0.85 – 2.59 0.174
Temperature [quadratic] 0.42 0.23 – 0.76 0.006 0.93 0.09 – 9.70 0.949 0.45 0.25 – 0.81 0.009
Isolate [C1] 1.58 0.59 – 4.28 0.366 1.49 0.55 – 4.03 0.432
Isolate [C18] 0.39 0.15 – 1.02 0.060 0.36 0.13 – 0.94 0.042
Isolate [C24] 0.70 0.26 – 1.83 0.464 0.62 0.23 – 1.66 0.346
DeathAge 1.11 1.07 – 1.14 <0.001 1.10 1.06 – 1.14 <0.001 1.09 1.06 – 1.12 <0.001
DeathAge × Temperature
[linear] 		0.96 0.92 – 1.01 0.107
DeathAge × Temperature
[quadratic] 		0.98 0.93 – 1.03 0.517
Observations 162 162 162
R2 Nagelkerke 0.348 0.362 0.299 
brant(modele14)
## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      16.1    18  0.59
## Temperature.L    4.04    3   0.26
## Temperature.Q    6.78    3   0.08
## IsolateC1    0.22    3   0.97
## IsolateC18   -461.22 3   1
## IsolateC24   2.63    3   0.45
## DeathAge 9.09    3   0.03
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds
# GSimplified stages
FC1 = polr(stade_simplifie ~ Isolate + Temperature * DeathAge, 
                data = infected, Hess = TRUE)
model_selection = dredge(FC1)

model_selection
## Global model call: polr(formula = stade_simplifie ~ Isolate + Temperature * DeathAge, 
##     data = infected, Hess = TRUE)
## ---
## Model selection table 
##    (Int)     DtA Isl Tmp DtA:Tmp df   logLik  AICc delta weight
## 8      + 0.09511   +   +          8 -125.533 268.0  0.00  0.435
## 16     + 0.08670   +   +       + 10 -123.565 268.6  0.58  0.325
## 6      + 0.08108       +          5 -130.010 270.4  2.40  0.131
## 14     + 0.07302       +       +  7 -128.273 271.3  3.27  0.085
## 4      + 0.08703   +              6 -131.054 274.7  6.64  0.016
## 2      + 0.07531                  3 -134.837 275.8  7.82  0.009
## 5      +               +          4 -148.833 305.9 37.91  0.000
## 1      +                          2 -151.702 307.5 39.47  0.000
## 7      +           +   +          7 -148.297 311.3 43.31  0.000
## 3      +           +              5 -151.140 312.7 44.66  0.000
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  stade_simplifie stade_simplifie
Predictors Odds Ratios CI p Odds Ratios CI p
Caul Grapes|Imm Mature 3.04 0.79 – 11.66 0.105 1.72 0.38 – 7.84 0.483
Imm Mature|2nd Caul 71.72 15.03 – 342.18 <0.001 43.91 8.21 – 234.87 <0.001
DeathAge 1.10 1.07 – 1.14 <0.001 1.09 1.06 – 1.13 <0.001
Isolate [C1] 1.54 0.57 – 4.15 0.392 1.43 0.53 – 3.87 0.477
Isolate [C18] 0.40 0.15 – 1.03 0.063 0.36 0.13 – 0.94 0.042
Isolate [C24] 0.77 0.29 – 2.05 0.604 0.69 0.25 – 1.84 0.458
Temperature [linear] 1.57 0.89 – 2.82 0.126 14.72 1.41 – 164.42 0.028
Temperature [quadratic] 0.45 0.25 – 0.82 0.010 1.51 0.15 – 15.81 0.728
DeathAge × Temperature
[linear] 		0.95 0.91 – 1.00 0.056
DeathAge × Temperature
[quadratic] 		0.97 0.92 – 1.02 0.318
Observations 162 162
R2 Nagelkerke 0.326 0.347 
best_mod = polr(stade_simplifie ~ Isolate + Temperature + DeathAge, 
                data = infected, Hess = TRUE)
tab_model(best_mod)
  stade_simplifie
Predictors Odds Ratios CI p
Caul Grapes|Imm Mature 3.04 0.79 – 11.66 0.105
Imm Mature|2nd Caul 71.72 15.03 – 342.18 <0.001
Isolate [C1] 1.54 0.57 – 4.15 0.392
Isolate [C18] 0.40 0.15 – 1.03 0.063
Isolate [C24] 0.77 0.29 – 2.05 0.604
Temperature [linear] 1.57 0.89 – 2.82 0.126
Temperature [quadratic] 0.45 0.25 – 0.82 0.010
DeathAge 1.10 1.07 – 1.14 <0.001
Observations 162
R2 Nagelkerke 0.326 
brant(best_mod) # Not much better, but on death age and omnibus, same concluisons, so I kept the other one.
## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      12.46   6   0.05
## IsolateC1    0.04    1   0.85
## IsolateC18   0.97    1   0.32
## IsolateC24   1.41    1   0.24
## Temperature.L    1.19    1   0.28
## Temperature.Q    0.72    1   0.4
## DeathAge 9.95    1   0
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds
tab_model(modele14)
  stade_max
Predictors Odds Ratios CI p
Cauliflower|Grapes 3.24 0.85 – 12.35 0.084
Grapes|Immature 3.60 0.94 – 13.75 0.061
Immature|Mature 4.40 1.14 – 16.94 0.032
Mature|2nd Cauliflower 95.60 19.81 – 461.41 <0.001
Temperature [linear] 1.53 0.87 – 2.73 0.145
Temperature [quadratic] 0.42 0.23 – 0.76 0.006
Isolate [C1] 1.58 0.59 – 4.28 0.366
Isolate [C18] 0.39 0.15 – 1.02 0.060
Isolate [C24] 0.70 0.26 – 1.83 0.464
DeathAge 1.11 1.07 – 1.14 <0.001
Observations 162
R2 Nagelkerke 0.348

Visualisation

## `summarise()` has grouped output by 'Temperature', 'Isolate'. You can override
## using the `.groups` argument.


The quadratic effect of temperature (OR = 0.42, p = 0.006) indicates a non-linear relationship, where parasite maturation is optimal at intermediate temperatures (30°C) and decreases at both lower (20°C) and higher (40°C) extremes. At extreme temperatures, the odds of reaching a more advanced maturation stage decrease by 58%. There is also a non significant trend of a less efficient maturation for C18, which is consistent with the previous analysis for the second stage cauliflower formation.

Predicted composition per temperature

6. Relationship between host castration and parasite fitness

This is a supplementary analysis that we decided not to include as it would be outside the scope of our article.

FC1 = glmmTMB(data=data_mat, mature_load/100000 ~ Never_repro * Temperature + Isolate ,family= gaussian(), REML=T) 

FC1R = glmmTMB(data=data_mat, mature_load/100000 ~ Never_repro * Temperature + Isolate + (1|Batch), family= gaussian(), REML=T) 

AICc(FC1, FC1R)
##      df     AICc
## FC1  10 1424.867
## FC1R 11       NA
data_mat$Nb_offspring = !is.na(data_mat$Nb_offspring)

FC1 = glmmTMB(data=data_mat, mature_load/100000 ~ Never_repro * Temperature * Isolate * DeathAge,
               family= gaussian(), na.action = na.fail)
model_selection = dredge(FC1)

head(model_selection)
## Global model call: glmmTMB(formula = mature_load/1e+05 ~ Never_repro * Temperature * 
##     Isolate * DeathAge, data = data_mat, family = gaussian(), 
##     na.action = na.fail, ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Nvr_rpr) cnd(DtA:Isl)
## 8     -48.740          +   1.9350        +       10.920             
## 40    -35.850          +   1.6910        +      -19.910             
## 56     15.510          +   0.6365        +      -28.080            +
## 4     -39.320          +   1.8480        +                          
## 24     -2.410          +   0.9804        +        9.321            +
## 20      9.883          +   0.8171        +                         +
##    cnd(DtA:Nvr_rpr) df   logLik   AICc delta weight
## 8                    7 -705.627 1426.1  0.00  0.303
## 40           0.6523  8 -704.766 1426.6  0.52  0.233
## 56           0.7908 11 -701.897 1427.8  1.71  0.129
## 4                    6 -707.590 1427.8  1.72  0.128
## 24                  10 -703.152 1427.9  1.87  0.119
## 20                   9 -704.605 1428.5  2.47  0.088
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  mature_load/1e+05 mature_load/1e+05 mature_load/1e+05 mature_load/1e+05 mature_load/1e+05
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
DeathAge 1.93 1.45 – 2.42 <0.001 1.69 1.08 – 2.30 <0.001 0.64 -0.48 – 1.75 0.264 1.85 1.36 – 2.34 <0.001 0.98 -0.06 – 2.02 0.065
Isolate [C1] 7.59 -7.98 – 23.15 0.340 6.74 -8.79 – 22.27 0.395 -60.72 -135.66 – 14.22 0.112 8.57 -7.18 – 24.32 0.286 -61.46 -137.05 – 14.12 0.111
Isolate [C18] 28.90 14.05 – 43.74 <0.001 28.18 13.38 – 42.98 <0.001 -46.60 -112.61 – 19.41 0.166 27.62 12.62 – 42.61 <0.001 -37.89 -103.59 – 27.81 0.258
Isolate [C24] 9.61 -5.91 – 25.13 0.225 9.10 -6.35 – 24.55 0.248 -31.03 -97.81 – 35.75 0.362 8.56 -7.14 – 24.26 0.286 -30.78 -98.14 – 36.58 0.370
Never repro 10.92 0.19 – 21.65 0.046 -19.91 -67.05 – 27.23 0.408 -28.08 -75.34 – 19.18 0.244 9.32 -1.34 – 19.98 0.087
DeathAge × Never repro 0.65 -0.32 – 1.62 0.188 0.79 -0.18 – 1.76 0.112
DeathAge × Isolate [C1] 1.46 -0.18 – 3.10 0.081 1.51 -0.15 – 3.16 0.074
DeathAge × Isolate [C18] 1.51 0.21 – 2.82 0.023 1.36 0.06 – 2.66 0.040
DeathAge × Isolate [C24] 0.84 -0.50 – 2.18 0.219 0.85 -0.50 – 2.20 0.219
Observations 145 145 145 145 145
R2 / R2 adjusted 0.395 / 0.368 0.402 / 0.371 0.425 / 0.382 0.378 / 0.356 0.415 / 0.376 
best_model = get.models(model_selection, 1)[[1]]

summary(best_model)
##  Family: gaussian  ( identity )
## Formula:          mature_load/1e+05 ~ DeathAge + Isolate + Never_repro + 1
## Data: data_mat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1425.3   1446.1   -705.6   1411.3      138 
## 
## 
## Dispersion estimate for gaussian family (sigma^2):  987 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -48.739     13.903  -3.506 0.000456 ***
## DeathAge       1.935      0.249   7.770 7.87e-15 ***
## IsolateC1      7.586      7.944   0.955 0.339599    
## IsolateC18    28.896      7.576   3.814 0.000137 ***
## IsolateC24     9.611      7.920   1.213 0.224975    
## Never_repro   10.923      5.476   1.995 0.046064 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.348 0.352 0.484 0.78 0.292 0.672 0.328 0.38 0.696 0.468 0.104 0.048 0.164 0.24 0.468 0.428 0.596 0.376 0.268 0.412 ...


It seems that there is a possible slight positive impact of castration on parasite load. However, the effective in the different categories are really small.

Visualisation

## `summarise()` has grouped output by 'Never_repro', 'Isolate'. You can override
## using the `.groups` argument.

## `summarise()` has grouped output by 'Never_repro', 'Isolate'. You can override
## using the `.groups` argument.

Session Info

sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Ventura 13.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Paris
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] brant_0.3-0         ordinal_2023.12-4.1 car_3.1-3          
##  [4] carData_3.0-5       MASS_7.3-64         brms_2.22.0        
##  [7] Rcpp_1.0.14         coxme_2.2-22        bdsmatrix_1.3-7    
## [10] glmmTMB_1.1.10      forcats_1.0.0       tidyr_1.3.1        
## [13] dplyr_1.1.4         gridExtra_2.3       survminer_0.5.0    
## [16] ggpubr_0.6.0        survival_3.8-3      MuMIn_1.48.4       
## [19] sjPlot_2.8.17       DHARMa_0.4.7        patchwork_1.3.0    
## [22] ggplot2_3.5.1       readr_2.1.5        
## 
## loaded via a namespace (and not attached):
##   [1] tensorA_0.36.2.1     jsonlite_1.8.9       datawizard_1.0.2    
##   [4] magrittr_2.0.3       estimability_1.5.1   farver_2.1.2        
##   [7] nloptr_2.1.1         rmarkdown_2.29       vctrs_0.6.5         
##  [10] minqa_1.2.8          effectsize_1.0.0     rstatix_0.7.2       
##  [13] htmltools_0.5.8.1    distributional_0.5.0 broom_1.0.7         
##  [16] Formula_1.2-5        sjmisc_2.8.10        sass_0.4.9          
##  [19] bslib_0.9.0          plyr_1.8.9           emmeans_1.10.7      
##  [22] zoo_1.8-13           cachem_1.1.0         TMB_1.9.17          
##  [25] commonmark_1.9.5     mime_0.12            iterators_1.0.14    
##  [28] lifecycle_1.0.4      pkgconfig_2.0.3      gap_1.6             
##  [31] sjlabelled_1.2.0     Matrix_1.7-2         R6_2.5.1            
##  [34] fastmap_1.2.0        shiny_1.10.0         rbibutils_2.3       
##  [37] digest_0.6.37        numDeriv_2016.8-1.1  colorspace_2.1-1    
##  [40] qgam_1.3.4           labeling_0.4.3       km.ci_0.5-6         
##  [43] abind_1.4-8          mgcv_1.9-1           compiler_4.4.2      
##  [46] doParallel_1.0.17    bit64_4.6.0-1        withr_3.0.2         
##  [49] backports_1.5.0      performance_0.13.0   ggsignif_0.6.4      
##  [52] sjstats_0.19.0       ucminf_1.2.2         loo_2.8.0           
##  [55] tools_4.4.2          httpuv_1.6.15        glue_1.8.0          
##  [58] promises_1.3.2       nlme_3.1-167         gridtext_0.1.5      
##  [61] grid_4.4.2           checkmate_2.3.2      generics_0.1.3      
##  [64] gtable_0.3.6         tzdb_0.4.0           KMsurv_0.1-5        
##  [67] data.table_1.16.4    hms_1.1.3            xml2_1.3.6          
##  [70] foreach_1.5.2        pillar_1.10.1        markdown_2.0        
##  [73] stringr_1.5.1        vroom_1.6.5          later_1.4.1         
##  [76] posterior_1.6.1      splines_4.4.2        ggtext_0.1.2        
##  [79] lattice_0.22-6       bit_4.5.0.1          tidyselect_1.2.1    
##  [82] knitr_1.49           reformulas_0.4.0     litedown_0.6        
##  [85] stats4_4.4.2         xfun_0.51            bridgesampling_1.1-2
##  [88] matrixStats_1.5.0    stringi_1.8.4        yaml_2.3.10         
##  [91] boot_1.3-31          codetools_0.2-20     evaluate_1.0.3      
##  [94] tibble_3.2.1         cli_3.6.3            RcppParallel_5.1.10 
##  [97] xtable_1.8-4         parameters_0.24.2    Rdpack_2.6.2        
## [100] munsell_0.5.1        jquerylib_0.1.4      survMisc_0.5.6      
## [103] ggeffects_2.2.1      coda_0.19-4.1        parallel_4.4.2      
## [106] rstantools_2.4.0     bayestestR_0.15.2    gap.datasets_0.0.6  
## [109] bayesplot_1.11.1     Brobdingnag_1.2-9    lme4_1.1-36         
## [112] mvtnorm_1.3-3        scales_1.3.0         insight_1.1.0       
## [115] purrr_1.0.4          crayon_1.5.3         rlang_1.1.5